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THE EXTENSION OF HOUSTON’S METHOD WITH 
APPLICATION TO DEBYE ®,' 


D. D. BEtTTs 


ABSTRACT 


Houston’s method of integration over the unit sphere of functions of cubic 
symmetry is extended by developing integration formulas using 9 and 15 values 
of the integrand. The 3-, 6-, 9-, and 15-term Houston formulas are used to 
calculate Debye 60 values for a number of cubic crystals. The convergence of 
the results is excellent except for the most anisotropic crystals. 


INTRODUCTION 


Houston (1948) proposed a method for integrating functions of complete 
cubic symmetry (corresponding to the identity representation of the point 
group O,) over the unit sphere using only three appropriately weighted values 
of the integrand. Further such formulas using up to 6 values of the integrand 
have been presented by Betts, Bhatia, and Wyman (1956) and used to com- 
pute Debye 9» values. Formulas for hexagonal, trigonal, and tetragonal 
symmetries have also been presented and used by Betts, Bhatia, and Horton 
(1956). Here we report on a further extension of Houston’s method for cubic 
crystals, motivated on the one hand by the “good”’ @o values obtained using 
the 6-term formula for most crystals and on the other hand by the still un- 
certain values obtained for the elastically very anisotropic alkali metals. 


DETERMINATION OF THE INTEGRATION FORMULAS 


Suppose [(0,¢) has complete cubic symmetry and we require 


2x 7 
(1) J= f d¢ | sin@d0 (6,4). 
70 0 


Houston’s method consists in expanding /(@,@) in the Kubic Harmonics of 
Von der Lage and Bethe (1947) having the same symmetry, which we denote 


by K in(9,) 
(2) I = >. A yn Kune 


l,m 


1Manuscript received October 3, 1960. ch 8 
Contribution from the Department of Physics, University of Alberta, Edmonton, Alta. 
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The degree of the polynomial is 2/ when expressed as is usual in terms of 
x, y, and gz, the direction cosines of 6 and ¢; m distinguishes K’s of the same 
degree. 

As all the other K;», are orthogonal over the sphere to Ko = 1 the integral 
(1) becomes J = 47A oo. 

In practice the series (2) uses only a small number of polynomials up to a 
certain degree, 2L. Then in (2) N special pairs of values of @ and @ (or N 
directions) are chosen for evaluating J and the K’s. The resulting NV equations 
determine the N unknowns Ajm. In particular A oo will be of the form 

N 


(3) Ago = 2» w,1(8in¢:). 

For example Houston’s original 3-term formula is 

(4) 35A 6? = 107(100)+16/(110)+97(111). 

The further formulas of Betts, Bhatia, and Wyman (1956) use values of the 

integrand in the 210, 211, and 221 directions, culminating in the formula 

(5) 1081,080 A‘? = 117,6037 (100) +76,5447 (110)+17,4967 (111) 
+381,2507 (210)+311,0407 (211)+177,147] (221). 


A natural way of choosing further directions in extending the method 
appears from examination of the plane projection of 1/48th of the unit sphere 
shown in Fig. 1. The procedure is to divide the large triangle into smaller and 





Fic. 1. Distribution over 1/48th of the unit sphere of points for evaluation of the integrand 
for use in Houston-type integration formulas. 


smaller triangles each time using all the corner points. Thus the first approxi- 
mation is formula (4), the third the 6-term formula (5) where the triangle is 
divided into four nearly similar triangles. Taking three further directions 
(411), (431), and (433) divides the original triangle into eight smaller triangles 
to give a 9-term approximation. Finally we obtain a 15-term approximation, 
by including the additional directions (410), (411), (480), (421), (432), and 
(443). The resulting formulas are 
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(6) A$? = 0.05637817I (100) -+0.04952710J (110) 
—0.05588614 (111)+0.17827058 (210) 
+0.102947717 (211)-+0.073098757 (221) 
+0.207116387 (411) +0.21376348I (431) 
+0.17478397 (433) 


(7) Ate’ = +0.04367057 (100) +0.0229221 J (410) 
+0.34002647 (210) —0.8133145J (430) 
+0.9623617 I (110) +0.2285255/ (411) 
— 0.5323867] (421)+2.2314783I (431) 
— 1.9823060/ (441) +0.8005275/ (411) 
—3.23705331 (432)+2.3075765/ (221) 
+2.5033260] (433) — 1.7484174I (443) 
—0.1269366J (111). 


The weights, w,, in (6) and (7) had to be found approximately by inversion 
of the matrix of K’s on the University of Alberta’s LGP-30 electronic digital 
computer.* The inverses were checked by multiplication of the calculated 
inverses by the original matrices, and the entries in the top row of the product 
matrix were found to be unity or zero respectively to at least six decimal 
figures. Thus the weights can be regarded as correct to seven figures. 


APPLICATIONS 


It is well known that Houston-type formulas can be used to calculate 00, 
the Debye characteristic temperature at 0° K. They have also been used by 
Horton and Schiff (1958) to calculate the coefficient B, in the low-temperature 
expansion of the Debye 9 given by 


(8) Or = Oo >. Bu(T/Oo)™ 


and by Horton (1960) to calculate an average Griineisen y in the theory of 
thermal expansion of solids. Although the tables of de Launay (1956) provide 
an accurate and quick way of computing Debye @o, they are not applicable 
to the calculation of any other physical quantities. 

*Actually the completely orthogonalized K’s need not be used but only polynomials of the 


form 
Bi = xttytmgte +x y2ng2l + x2ny2lz2m +-x2lytngim 2m y2izen + x2ny2mz2! ‘and bimn 


where Dima is chosen to make the integral of Rim, over the unit sphere vanish. The number 
of k’s of a given degree to be used is determined by formula (A.10) of Betts, Bhatia, and 
Wyman (1956). The evaluation of the 6’s is trivial. The latter polynomials were the ones 
used to calculate formula (7). 
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We have used here formulas (4)—(6) to calculate successive approximations 
e®, 6, and 6%) for 22 crystalline solids from the well-known formula of 


Debye 


_4(9 _» ¥- 
(9) ws (2 WAca 


where h and & are Planck’s and Boltzmann’s constants, p is the density, and 
M the average mass of an atom of the substance. Ago is the average over 
the unit sphere of 





(10) (6, $) = x {v4(0, )}~° 


where v, are the elastic wave velocities. For seven of these substances in- 
cluding the most anisotropic we have also calculated the 15-term approxima- 
tion, 0“. Table I contains the results. 

As a further check on the formulas (4)—(7) we have calculated approximately 


the integral 
2r a : 
(11) A= J do | sin 6d (arc sin*x+-arc sin*y+-arc sin’z) 
0 0 
= 3x(n°—8). 
The results are presented in Table II. 
TABLE II 


Approximations to A = f2*d¢fj sin 6 déx 
[arc sin? x+-arc sin? y+-arc sin? zs] 








1 
(number of values of 





integrand used) Aw 
3 19.619 
6 18.096 
9 17.797 
15 17.739 


Exact 17.621 


DISCUSSION AND CONCLUSIONS 


For substances with anisotropies from .5 to 2.0 the 6° and 6° approxi- 
mations to Qo agree to closer than 1/10th of a per cent. Hence we did not, 
in general, bother for such substances to calculate 6“. For germanium, how- 
ever, with 7 = 1.68, 6) and 6“ agree to one part in the fifth figure. Even 
for such an anisotropic substance as lead 0‘ and 6“ agree to better than 
two one-hundredths of a per cent. It is somewhat disappointing to see then 
that for sodium and lithium the 9-term approximation is no closer to the 
15-term approximation than to the 6-term one. Probably a tremendously 
large number of terms in the harmonic expansion of J(6,¢) in these cases 
are needed to represent the function accurately. 








238 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


Table II shows the satisfactory convergence of Houston’s method for 
integration of a function whose series expansion is known to converge very 
slowly and for which the exact value of the integral is known. 

We conclude that the present extension of Houston’s method provides 
sufficiently accurate integration formulas for computing Debye Qo of all 
substances but alkali metals. Although it would not be difficult to extend the 
method systematically to still higher degrees of approximation it is doubtful 
if the labor is worth while. The formulas may of course be used for the precise 
integration over the unit sphere of any function of complete cubic symmetry. 
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ASYMMETRY IN THE RECOVERY FROM A VERY DEEP 
FORBUSH-TYPE DECREASE IN COSMIC-RAY INTENSITY! 


D. C. Rose AnpD S. M. LAPOINTE? 


ABSTRACT 


The intensity-time curves for cosmic rays recorded at some 30 stations dis- 
tributed all over the world are examined for structure in the recovery period from 
the third in a series of three closely spaced Forbush-type decreases which occurred 
in the middle of July 1959. It is shown that the structure of intensity peaks is 
regular and that these occur at each station at the same effective local time. 
It is found that this is consistent with the hypothesis that recovery from a very 
deep Forbush-type decrease is first apparent in directions making 15° and 165° 
with the sun-earth line respectively. The analyses suggest further, that during 
recovery from this deep Forbush decrease temporary openings appeared in the 
intensity depressing mechanism which allowed intensity increases in limited 
directions. 


t I. INTRODUCTION 

In the middle of July 1959 there occurred in close succession three large 
Forbush-type decreases. The complete intensity curves for that period at the 
six stations—Ottawa, Sulphur Mountain, Churchill, and Resolute, operated 
by the National Research Council of Canada; Deep River, operated by Atomic 
Energy of Canada Limited; and Thule, operated in close co-operation with 
the Canadian stations by the Bartol Research Foundation—have already been 
published (Carmichael and Steljes 1959; Wilson, Rose, and Pomerantz 1960). 
The first large Forbush decrease occurred on July 11 followed by a very slow 
recovery. Before recovery to the pre-event level a second decrease occurred 
on July 15. This was followed by an unusually rapid increase which has been 
interpreted by Steljes and Carmichael (1960) as a beam of solar injected 
particles. After this event, late on July 17, there occurred a third Forbush 
decrease bringing the intensity down toa very low point indeed. The recovery 
from this last decrease reveals some interesting structure. Among other features, 
attention has already been called to a peak in the hourly values of the intensity 
early on July 18 and marked D in the curves shown by Wilson, Rose, and 
Pomerantz (1960) and similarly marked in Fig. 1 of this paper. We would 
like here to attempt a detailed analysis of this peak and following events 
(E and F, Fig. 1) making use of data from a large number of stations.* This 
analysis will be limited mostly to variations in the nucleonic component of 
cosmic rays as measured by standard neutron monitors of the type designed 
for the International Geophysical Year. 


1Manuscript received November 15, 1960. 

Contribution from the Division of Pure Physics, National Research Council, Ottawa, 
Canada. 
Issued as N.R.C. No. 6133. 

2Department of Physics, University of Montreal, Montreal, Quebec. 

*The data from about 30 stations in all were used. These data were obtained in most cases 
directly from the laboratories concerned by data exchange arrangements during the I.G.Y. and 
1.G.C. Carmichael and Steljes made a world-wide collection of data for this event and since 
our analysis of this event was carried out in close co-operation with them, in a few cases where 
we did not have the data ourselves their collection was used. In the final paragraph of this 
paper we have expressed appreciation to our many colleagues for making these data available. 
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Fic. 1. Intensity-time curves on July 14 through 20, 1959, at: Ottawa, Churchill, Sulphur 
Mountain, College, and Mt. Wellington. Increasing geographic longitude to the west reading 
down. Average line through flat portion before the decrease on July 15 arbitrarily taken as 
100%; other horizontal lines represent 5% variations. The time is universal time. The three 
peaks concerned here are marked D, E, and F, on the Churchill curve. They are clearly apparent 
on the other curves, 
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2. TIME OF OCCURRENCE OF PEAK “D” 


The first point of interest is to determine where and when this sharp event 
occurred. Figure 1* shows the intensity-time curves from July 14 to July 20 
at some stations at which the event is clearly present; they are arranged reading 
downward in order of geographic longitude increasing to the west. The time is 
universal time. An average line through the level portion before the decrease 
on July 15 arbitrarily represents 100% intensity; other horizontal lines 
represent 5% variations. For comparison Fig. 2 shows the intensity-time 
curves at other stations at which no sharp event is apparent early on July 18. 
These are divided in two groups; group A the high and low latitude stations, 
and group B the intermediate latitude ones. The stations in group B are also 
arranged in order of increasing geographic longitude to the west reading down. 
Figure 3 shows a map of the locations at which the early event ‘‘D” (Fig. 1) on 
July 18 did occur and those at which it did not. The positive events are seen 
to be well-grouped inside the oblong curve. 

Let us consider next, the time at which this event occurred at the positive 
stations. The characteristic time is chosen as that represented by an intensity 
weighted mean time. This, expressed in universal time, as well as station local 
time, is shown in Table I together with the stations’ geographic co-ordinates. 


TABLE I 
Intensity weighted mean time of occurrence of the sharp early event 
“D” on July 18, 1959, at 12 intermediate latitude stations (10 north, 
2 south) expressed in U.T. (universal time) and L.T. (station local time). 
Increasing longitude to the west reading down 











Geographic 

Station Lat. Long. Tu.T. TL. Tt: 
Ottawa 45° N. 76° W. 0415+0005 23.3 
Deep River 46° N. 78° W. 0505 +0005 23.9 
Chicago 42° N. 88° W. 0530 +0030 23.6 
Churchill 59° N. 94° W. 0555 +0005 23.6 
Lincoln 41° N. 97° W. 0555 +0100 22.5 
Climax 39° N. 106° W. 0615-+0030 23.2 
Sulphur Mt. 51° N. 116° W. 0645 +0005 23.1 
Berkeley 38° N. 122° W. 0600 +0100 21.9 
College 65° N. 148° W. 0830 +0030 22.6 
Sydney 34° S. 151° E. 0915+0100 19.4 
Mt. Wellington 43° S. 147° E, 1045+0030 20.6 
5 


Mt. Norikura 36° N. 138° E. 0915+0100 18. 








It might be noted that the time at the five Canadian stations is given with less 
uncertainty than elsewhere since it was determined from 10-minute values of 
the intensity. The uncertainty is arbitrarily given as one-half the smallest 
interval of the intensity plot used. This amounts to 5 minutes for the Canadian 
stations, 30 minutes for stations giving hourly values, and 1 hour for stations 
giving bihourly values. The most striking feature of this table is the gradual 
increase of the characteristic time of the event (when expressed in U.T.) with 


*We are grateful to Carmichael and Steljes for the original plotting of these curves. They 
have developed automatic plotting and normalization techniques which are very useful. 








242 CANADIAN JOURNAL OF PHYSICS, VOL. 39, 1961 


olf { mene 


90 | + +f 
| 
4 5 
im! | ul 
100 dt Al, %/ 
a k Ay HAWAII 
t ||. fh Nga * 
MN 95 mt Ml * | fea 
— 
a 1 
5 1 hol 
— r ye © 4 hy 
ey \! KAMPALA 
= v A n 
w han ) 
2 3 ae 
ml 
iw | 
, (y 
uw loony MINA 
o i : yi] AGUILAR 
> \ } 
p\ Vv 


14 OIE 6 , \8 9 2 
DAYS UT. 
JULY 1959 


4 15 16 I7 18 19 20 


} UPPSALA 
95 \ fh Ls 
‘| 


90 


LI f ri 


vr 
100 8 + -y HERSTMONCEUX 
i fo 
95 } J\ nf | 
| ‘ay | | 
1a 
90! +4 } j 
Y I 
ll 
loom A 4 ELLSWORTH 
95 } { A 
} 4 (rat 
| Ay 
iy | i) 
Wr VAN 
10€ I ‘} f 
1oorwy | % l MOUNT 
J WASHINGTON 
35 tt 
) nl 
; ft 


5 6 (7 18 19 20 


DAYS UT. 
JULY 1959 


Fic. 2. Intensity-time curves on July 14 through 20, 1959. Group A, at Thule for high 
latitudes, Hawaii, Kampala, and Mina Aguilar for low latitudes. Group B, at Uppsala, Zug- 
spitze, Herstmonceux, Ellsworth, and Mt. Washington for intermediate latitudes. Same 
co-ordinates as Fig. 1. 
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LOCATION OF THE EVENT EARLY ON JULY 18,1959 
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Fic. 3. Map of the stations at which the sharp event ‘“D” did and did not occur early 
on July 18, 1959. X or » means the event did occur in a neutron monitor or meson telescope 
respectiv: ely. O or » with a stroke through it means the event did not occur. The oblong curve 
encircles the main region where the event occurred. The abbreviations for the station names 
are shown in Tables II and III. 


increasing longitudes to the west, and the almost constant value of the mean 
time when expressed in local time. This indicates strongly that the beam of 
particles giving rise to the event did not arrive simultaneously from all direc- 
tions on earth, but rather that it was limited in size, fixed in space, and that it 
appeared at a given station as the station passed through it due to the earth’s 
daily rotation. 

The fact that stations limited to the region shown by the oblong curve in 
Fig. 3 are the only ones showing this peak also indicates that primarily the 
beam was only of a few hours’ duration. More details will be presented in 


Section 7 below. 


3. INTEGRAL INTENSITY SPECTRUM IN PEAK “D” 

The total intensity of the event measured at each one of the positive stations 
has been determined in the following way. First the intensity at each station 
has been expressed in per cent of the average of the intensity during the first 
9 hours of July 15. Then a baseline was arbitrarily chosen in each case as that 
line which represented the smoothed intensity without the peak. The intensity 
in each time interval was then taken as the difference between the actual 
intensity and this imaginary baseline; the resulting intensity is a figure which 
was less than 10% for each time interval. The sum of these intensities for each 
10-minute interval over the whole duration of the peak is taken to represent 
the total intensity. This total intensity is given in Table II for each one of the 
positive stations, together with the conventional geomagnetic co-ordinates and 
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TABLE II 


Total intensity of the early July 18, 1959, event ‘‘D” at 11 stations; the unit is explained 

in the text. The first two and the last column give the conventional geomagnetic co- 

ordinates and the conventional cutoff energy of vertically incident protons at each 
station. Increasing geomagnetic latitude reading down 


Abbr. Geomagnetic (degrees) Cutoff energy 
used in ————________—_- for protons in 
Station Figs. Lat. Long. Total Bev 
Mt. Norikura Mt. No 26 204 54 8.9 
Sydney Sy —42 226 101 3.64 
Berkeley Be 44 298 83 3.12 
Climax Cl 48 315 130 2.14 
Mt. Wellington Mt. We 52 225 148 1.36 
Chicago Chi 53 337 82 1.21 
Ottawa O 57 351 151 0.65 
Deep River DR 58 349 122 0.54 
Sulphur Mt. Su Mt. 58 300 252 0.54 
College Co 65 257 81 0.12 
Churchill Chu 69 323 195 0.04 





conventional cutoff energy for vertically incident protons at the station in 
question (centered dipole field). Although the event is not simultaneous at all 
stations, and although the beam intensity may vary with universal time and 
therefore also from station to station as they pass under the beam, nevertheless 
a plot of the intensity versus vertical cutoff energy may give some indication 
of the nature of the energy spectrum of the incoming particles. This has been 
done in Fig. 4. There is considerable variation in intensity with universal time 
but a comparison of stations with almost the same geographical longitude 
should give some indication of the shape of the spectrum when compared with 
a curve plotted from measurements on the ordinary total cosmic-ray spectrum. 
Such pairs of stations are Chicago and Churchill, Sulphur Mountain and 
Berkeley, Sydney and Mt. Wellington. Examine these, particularly in relation 
to the ordinary cosmic-ray spectrum represented by the solid line and, in fact, 
the general distribution of points suggests that the spectral distribution of 
peak ‘‘D” is somewhat flatter or at least no steeper than that for the ordinary 
cosmic-ray spectrum. The curve for the ordinary cosmic-ray spectrum was 
derived from Dorman (1957). Arbitrary units are used, therefore, the com- 
parison is only relative. The consistent high intensity of this peak being limited 
to the stations in the Western Hemisphere (enclosed by the closed curve in 
Fig. 3), as each station passed through the zone of maximum activity, suggests 
that this beam or flow of particles lasted for only about 10 hours. 

In view of the above evidence, we conclude that the beam giving rise to 
sharp peaks early on July 18, 1959, was made up of ordinary cosmic rays 
leaking through a “‘hole’’ in the modulating mechanism producing the Forbush 
decrease. For reasons mentioned in the last paragraph, comparisons of inten- 
sities in this way are subject to some criticism but there is additional evidence 
that the above conclusion is correct. An examination of the meson component 
shows that the peak ‘“‘D” appears in this component at an intensity consistent 
with this conclusion. The times of occurrence at a number of stations where 
meson data have been studied are the same as that for the nucleon component 
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INTEGRAL SPECTRUM 
OF THE EVENT EARLY ON JULY 18,1959 





















aera —- 
a | ORDINARY COSMIC RAY 
= SPECTRUM 
Zz 300}+——— 
= 
7 ° Chu 
<q 
c 
= 
— 
= 100o-————— 
— *Co 
> 
= 
wn 
= 
W 
5 
= 30 
ad 
a 
b 
°o 
e 
[ieee anenrmesl aegis Resi - ae 
Ol 0-3 | 3 10 


CONVENTIONAL VERTICAL CUTOFF ENERGY 
FOR PROTONS E, IN BEV 


Fic. 4. Total intensity of the early July 18, 1959, event ‘‘D” at various stations as a 
function of conventional cutoff energy for vertically incident protons at that station. The solid 
line represents the integral intensity spectrum (arbitrary ordinate) of ordinary cosmic-ray 
pratense (Dorman 1957). The abbreviations used for station names are shown in Tables II 
an , 


and the intensity, compared to the intensity of the nucleon peaks, is about that 
expected if the spectrum is that of the usual galactic cosmic rays. Ordinarily, 
the depth of a Forbush decrease is appreciably greater in the case of the 
nuclear component than for the meson component. The mean ratio for high 
latitude stations is about 2 though there is considerable variation (Fenton, 
Fenton, and Rose 1958). In the present case the depth of the Forbush decrease 
on July 17 is complicated by the large increase in the flux of soft particles 
occurring for a few hours before the July 17 decrease. The Forbush decrease on 
July 15, however, seems normal. The ratios of the depths of the decrease in 
nucleon intensity to the meson decrease for this event were 2.5, 2.3, and 2.1 
respectively for Ottawa, Sulphur Mountain, and Churchill. The corresponding 
ratios for the peak ‘“‘D”’ at these three stations are 1.0, 1.8, and 3.4 (obtained 
by numerical integration as explained above). These are within the ranges of 
expected values if the incoming particles have, at least approximately, the 
spectrum of normal galactic cosmic rays. 


4. DIRECTIONAL SENSITIVITY OF THE DETECTORS 


Because of the deflection of cosmic-ray particles in the earth’s magnetic field 
their apparent impact direction at the earth’s surface is not in general the 
direction in which they were travelling before entering the earth’s magnetic 
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field. In a given impact direction particles of different energies come from 
different directions often spoken of as the asymptotic directions. A neutron 
monitor for instance detects particles arriving within a cone of half angle 
about 50° centered on the vertical, with sensitivity decreasing rapidly with 
increasing angle from the vertical. Taking all these factors into account, one 
can determine the average asymptotic direction from which particles of a given 
energy spectrum reach the recording equipment. This has been called the 
effective direction (Dorman 1957) or the direction of maximum response 
(Fenton, McCracken, Rose, and Wilson 1959). Most particles come from a 
direction east of the station (Brunberg and Dattner 1953). The effective 
direction is also east of the station, the amount depending on the station’s 
geomagnetic latitude. We have calculated the effective directions of a neutron 
monitor covering all geomagnetic latitudes from the equator to the pole. The 
method of calculation and the results will be published elsewhere. From these 
calculations the effective directions are available for all the stations considered 
here and are listed in Table III. The accuracy of the calculations is adequate 


TABLE III 


Direction of maximum response to the east of the observing meridian expressed 
both in degrees and in time 


Geomagnetic Shift to Time 





Abbr. used lat. the east increment 

Station in Figs. (degrees) (degrees) (hr) 
Mt. Norikura Mt. No. 26 71 4.7 
Sydney Sy —42 79 5.3 
Berkeley Be 44 62 4.1 
Climax Cl 48 56 3.7 
Lincoln Li 51 47 3.1 
Mt. Wellington Mt. We —52 58 3.9 
Chicago Chi 53 43 2.9 
Ottawa O 57 36 2.4 
Deep River DR 58 38 2.5 
Sulphur Mt. Su Mt. 58 28 1.9 
College Co 65 16 Det 
Churchill Chu 69 19 1.3 
Mawson Ma —73 7 0.5 
Uppsala Up 59 54 3.6 
Munich Mu 49 70 4.7 
Zugspitze Zu 48 74 4.9 
Pic du Midi PM 46 78 5.2 
Herstmonceux He 54 60 4.0 
Leeds Le 57 55 3.7 
Ellsworth El —67 24 1.6 
Rio de Janeiro Rio —13 78 5.2 
Mina Aguilar M Ag —12 82 5.5 
Mt. Washington Mt. Wa 56 39 2.6 
Resolute Re 83 —1 —0.1 





for the purposes considered here. In Table III the directions are expressed as 
a longitude shift to the east of the observers’ meridian, that is, the station 
passes under the effective direction at a time later than the time of observation 
of the maximum by the amount shown in Table II]. The corrected times or the 
times the stations are in the meridian of the effective direction are shown in 
Table IV. These are shown in local time and within the order of over-all 
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TABLE IV 


Intensity weighted mean time of event ‘D” 

early on July 18, 1959, in the nucleonic com- 

ponent, corrected for direction of maximum 
response, and expressed in local time 





TL. T1.1. 

Station actual corrected 
Ottawa 23.3 Ba 
Deep River 23.9 2.4 
Chicago 23.6 2.5 
Churchill 23.6 0.9 
Lincoln 22.5 1.6 
Climax 23.2 2.9 
Sulphur Mt. 23.1 1.0 
Berkeley 21.9 2.0 
College 22.6 23.7 
Sydney 19.4 0.7 
Mt. Wellington 20.6 0.5 
Mt. Norikura 18.5 23.2 
Average 1.3 


accuracy; the local time is almost the same at all stations. The average is 1.3 
hours. The four stations from which most reliable times can be derived namely 
Ottawa, Deep River, Churchill, and Sulphur Mountain, give an average 
corrected local time of 1.5 hours with extremes of +0.9 and —0.6 hours. This 
corresponds to a direction in the sky making an angle of about 165° with the 
sun-earth line. The opening in the modulating mechanism which in the case 
of peak ‘‘D”’ lasted about 10 hours occurred on the side of the earth almost 
opposite the sun. 


5. ASYMMETRY IN THE RECOVERY FROM THIS DEEP FORBUSH EVENT 


Although this hole in the modulating mechanism only appeared open for 
about 10 hours when it was first apparent, there is evidence that recovery in 
this direction remained more favorable (allowing varying intensities of passage 
of recovery particles) for at least two rotations of the earth. Further, another 
direction of favored recovery appeared in a direction almost opposite the 
above or nearly in the sun-earth line. The time-intensity curves in Figs. 1 and 
2, in many cases reveal more peaks later on in the day on July 18 and on July 19 
during the early part of the recovery from the Forbush decrease. These peaks 
are listed in Table V, when they are clearly distinguishable, together with the 
approximate time of their occurrence expressed in universal time: the un- 
certainty in this time, read directly off the intensity-time graphs is taken as 
+2 hours. These events are then plotted in Fig. 5 after correction for direction 
of maximum response. This graph shows geographic longitude at which the 
event occurred as a function of corrected universal time when it occurred. 
In other words, a point represents a direction in space under which the station 
concerned would be at the indicated universal time when the event occurred 
for that observer. The sloping lines show constant local time or constant 
direction with respect to the sun-earth line. It is remarkable how the events 
line up on slopes of constant local time situated 12 hours apart. This means, 





248 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


TABLE V 


List of the events on July 18 and 19, 1959. The time of occurrence (uncorrected) 

is expressed in U.T. A dash means no event is visible on the record. A question 

mark means data are incomplete in that region. The columns represent constant 

local time and, therefore, constant direction in the sky. Increasing longitude to 
the west reading down 





Station BLT. SObEA. 13.L.f. 22 L.T. 13 L.T. 


Ottawa 0415 14 24 16 
Deep River 0505 — 24 — 
Chicago 0530 16 24 —_ 
Churchill 0555 18 2 18 
Lincoln 0500 16 2 = 
Climax 0615 _ 1 ? 
Sulphur Mt. 0645 17 4 — 
Berkeley 0600 15 3 — 
College 0830 21 9 21 
Sydney 0915 . 6 ? 
Mt. Wellington 1045 23 9 — 
Mt. Norikura 0915 24 — 
Mawson 15 22 10 20 
Uppsala —_ —_ 8 20 
Munich — — 7 19 
Zugspitze — — 8 19 
Pic du Midi 13 -- — 17 
Herstmonceux — _— 7 20 
Leeds — _ 9 22 
Ellsworth 14 22 13 23 
Rio de Janeiro 12 -— a _— 
Mina Aguilar 13 — — ~- 
Mt. Washington — 24 14 1 





according to our hypotheses, that there existed two holes in the modulating 
mechanism diametrically opposed or 180° apart, one at 0100 L.T. and the 
other at 1300 L.T., making respectively 165° and 15° with the sun-earth line. 
These lasted for almost two complete revolutions of the earth and disappeared 
only on July 20 at which time the recovery from this last Forbush decrease 
seemed reasonably complete. However, since all the peaks do not occur at all 
stations as they pass under the direction of the opening the holes were not 
continuously open to the maximum extent. The first peak with an average 
position in time of 01.3 hours appeared to be fully present for only about 10 
hours represented by the width (in longitude) of the oval in Fig. 3. The length 
of this period is open to a question of interpretation but it is very clear that 
this peak was either missing or quite smali or possibly displaced in time as 
shown by the gap in points on or near the 01-hour line in Fig. 5 in the longitude 
region from 60° east to 70° west. The second and third peaks on the lines 
representing 13 and 22 hours local time seemed more uniformly distributed 
around the earth but were, of course, wider and in some cases sufficiently 
diffuse to make interpretation not entirely definite. 


6. COMPARISON WITH OTHER ASYMMETRIES IN THE FORBUSH DECREASE 
MECHANISM 


It has been reported by Fenton, McCracken, Rose, and Wilson (1959) that 
during a Forbush-type decrease, the intensity is first depressed in directions 
running from 30° to 120° to the west of the sun—earth line. The median direction 


aS AEE ATION 52 OED THR RE Tat 





MS 











ROSE AND LAPOINTE: COSMIC-RAY INTENSITY _ 249 























an a Poy / | ron / reed / fy 
Ww 70 ——— i Or i MW 
a ~ --— = M. Ag 
sof / / pee 
ee eet ee SSS =e? 
30;-— 
WI0}— - 
pn eae — Se coenee ean —— — bai os = He, PM. 
YO le the 
On ee ee ae 2 eee _o Up 
30}— / 
50}— / 
a ee a + Ma 
70}- 
‘i 
os Kj : 
AK i ~ ie 
Ww = 90, as o 
a D ~ ae 
eS nv & £& > 
FE g 5 Ss g 
ss. ™ 5 = = 
2 z = fe) Ny 
oO 1 So ~ SY] 
aa | s0~- = 
corniea ons — el — Mt.No 
2 - - at paki - We 
a 1SO}-— = eee - ee — Sy 
ra © 
S E 170}— vl 
WI70r So 
= 
2 
150/— He —_s aa #4 —Co 
130}— / / / / 
—— —_—— ——— Be 
~ _—s —— - Su Mt. 
HOF / / / / 
> ; cee ] Cli 
li 
—— — . ey = = ——i Chu 
90);— ——__" yy - | — - Chi 
/ / / / OR 
SF | > re Gg 
fa? ole bee Secon Sel | | beh] 
0 + 8 4 16 20 24 o 8 12 16 20 24 HOURS U.T. 
JULY 18,1959 JULY 19,1959 


Fic. 5. Time-occurrence diagram of the events recorded on July 18 and 19, 1959, at a 
number of stations. Ordinate: geographic longitude increasing to the west upwards; abscissa: 
universal time. Sloping lines represent constant local time, and hence constant direction in 
space. The abbreviations used for station names are shown in Tables II and III. 


to be first depressed, therefore, makes an angle of 75° with the sun-earth line 
and hence is perpendicular to the direction in which the intensity is first to 
recover, as determined above. The same authors give the median direction 
in which the intensity is last depressed as the one making an angle of 45° to 
the east of the sun-earth line; this by the way is also almost the direction of 
the maximum of the average daily variation at Hobart during 1957. All these 
directions are illustrated in Fig. 6. It might be noted that the directions of 
first and last depression in the Forbush mechanism were determined in the 
paper just quoted, from the analysis of 22 events and therefore seem to be a 


general enough feature. 
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Fic. 6. Characteristic directions with respect to the sun-earth line in the Forbush-type 
decrease mechanism. 


It seems very natural in a general sort of way, there being no definite model 
for the Forbush-type decrease mechanism, that the intensity should recover 
first in directions perpendicular to that in which it was first depressed. We 
therefore conclude from the foregoing analysis: 

(a) The early event on July 18, 1959, was not simultaneous at all the stations 
observing it on earth. 

(6) The energy spectrum of the incoming particles strongly suggests that 
they consist of ordinary cosmic rays. 

(c) This being the case, particles making 15° and 165° with the sun-earth 
line are first to leak through the depressing mechanism. 


7. DISCUSSION 


The night side (01 h L.T.) “hole’’ opened up first just over Ottawa. This 
onset seems to be very sharp since nearby Mt. Washington and M.I.T. 
(mesons) had just passed and did not record it at least not with any appreciable 
strength. Mt. Washington picks it up the next time around. 

The most outstanding and really the only significant gap is the one for the 
European stations around 24 U.T. July 18. It would seem that the night side 
hole was closed when they moved into position. Unfortunately, we have no 
data for the stations in Asia between the Australian and the European stations 
to follow the development of this intensity variation with time. But the 
openings certainly exist over the European stations in the next two positions. 

Some further features of these intensity increases which have not been 
emphasized yet merit discussion. The beam is assumed to be composed of 
ordinary cosmic-ray particles, but it is not detected at low latitude stations. 
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Our method of deciding whether or not the event occurred is necessarily coarse; 
and because of intensity fluctuations only relatively large cases can be dis- 
tinguished. Consider Sydney almost the lowest positive station, where the 
geomagnetic latitude is —42° and the vertical cutoff energy for protons is 
3.6 Bev. The next positive station available is Mt. Norikura where \ = 25° 
and E, = 8.9 Bev. At Sydney the maximum intensity in the peak is 4.5%. At 
Mt. Norikura it would be 4.5% X(3.6/8.9)'" or about 1.5%; this is barely 
detectable as an event. Therefore, the absence of events at lower latitude 
stations is consistent with an ordinary cosmic-ray beam leaking through holes 
and resulting in the events analyzed here. 

Two polar cap stations have been considered, namely, Thule and Resolute. 
One might expect the intensity peaks there to be as big as they are near the 
positive stations of hizhest latitude. However, the holes are probably limited 
in latitude, up to this point in the analysis only longitudinal asymmetries have 
been considered. 

According to our calculation of effective directions for neutron monitors, the 
asymptotic latitude of maximum response at Resolute is 64° N. geographical, 
and at Thule lies between 74° and 84° N. geographical. Hence, the absence of 
signal at the polar cap stations is consistent with the simple assumption that 
the holes extended to no more than 60° in the active meridian plane. If one 
takes 6 hours as the typical duration of the event (width of the beam), this 
means a divergence of 90° in the equatorial plane. If one were to assume the 
same size or divergence in a meridian plane, then one would conclude that the 
beam does not extend to more than 45° N. and 45° S. in asymptotic latitude. 
Therefore, the stations Thule and Resolute may have had no connecting orbits 


with the holes. 
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THE EMISSION SPECTRUM OF THE CCl RADICAL! 


R. D. Gorpon? AND G. W. KING 


ABSTRACT 


A rotational analysis of the 2780 A emission band obtained in a microwave 
discharge through CCl, vapor and photographed on a 20-ft grating spectrograph 
shows that a 2A;(b) — II,(a@) transition of the CCI radical is responsible, not 
2 — %[I(a) as reported by previous workers. Molecular constants are given for 
the combining states, as well as a vibrational analysis that identifies the 2780 A 
band as the (0-0) band. 


INTRODUCTION 


A high-frequency electrical discharge through carbon tetrachloride vapor at 
low pressure produces a series of emission bands in the 2710-2940 A region. 
Under low resolution four bands, lying at approximately 2714, 2780, 2850, 
and 2920 A, are prominent, those at 2780 and 2850 A being double-headed 
and especially intense in emission (Fig. 1). Various vibrational analyses have 
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Fic. 1. Microdensitometer profiles of CCl bands at (a) 2714, (b) 2780, (c) 2850, (d) 2920 A 
under low resolution. 
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been attempted (Asundi and Karim 1937; Horie 1939), the most extensive 
being that of Venkateswarlu (1950) on the basis of a ?2 — 7II electronic 
transition of the CCl radical. Kuzyakov and Tatevsky (1959a) agreed that 
this transition was responsible for the bands, but criticized inconsistencies in 
Venkateswarlu’s analysis, which gave a high value of 22 ev for the lower 
state dissociation energy. They proposed a new vibrational scheme, and also 
(19595) gave rotational analyses of two branches in each of the stronger 
bands of the system that supported a ?2 — *II(a) electronic assignment. How- 
ever, no vibrational analysis has been free of inconsistencies so far, mainly 
owing to uncertainty in identifying the observed band heads. 

We have rephotographed the spectrum under much higher resolution than 
was used in previous work, and a detailed study of the band at 2780 A shows 
that the transition is, in fact, ?A(6) — *II(@). This analysis enabled positive 
identification of heads and prominent features in the other bands to be made, 
resulting in a satisfactory vibrational scheme which shows the 2780 A band 
is the (0-0) transition. Both the rotational and vibrational analyses are con- 
sistent with CCI being the diatomic species responsible for the spectrum. 


EXPERIMENTAL 


Spectra were excited using a Raytheon ‘“Microtherm” model CMD 4, 
delivering 100 watts of microwave radiation at 2450 Mc/s. The antenna 
was mounted 1 cm from the outside of a 12-mm hard-glass discharge tube 
which was viewed end on by the spectrograph through a silica window. Helium 
gas at 0.1 mm pressure was flowed through the tube as a carrier gas, together 
with carbon tetrachloride at sufficient vapor pressure to change the red 
helium discharge to a bluish-white color; the carbon tetrachloride was of 
reagent grade quality, resublimed three times under vacuum. The visible 
discharge was some 5cm long. Low-resolution spectrograms were taken in 
the second order of a 1.5m Bausch and Lomb spectrograph model 11. High- 
resolution spectra were obtained on a 20-ft Ebert grating spectrograph (King 
1958), used in the second order with a resolving power in excess of 200,000. 
Using Kodak 103a0 plates, exposure times were about 4 hours, but increased 
with the age of the discharge tube, whose walls became coated with brown 
polymer. Plates were measured in both directions on a travelling microscope, 
each line being measured three times during each run. Iron arc spectra were 
used for calibration, with wavelengths taken from the M.I.T. tables (Harrison 
1939), together with the vacuum corrections of Edlén (1953). Vacuum wave- 
numbers and assignments for all the rotational lines of the 2780 A band are 
given in Table I. 

When the spectrum of pure helium was excited under identical conditions, 
plates of 11-hours’ exposure showed a very faint semicontinuous background, 
with no strong atomic lines in the region of the CCl spectrum. The spectra 
were free of contamination from this source, but the weak bands at 2714 and 
2920 A, which were only photographed under low resolution, were, apart 
from heads, just slightly more intense than a continuous background ascribable 
to the chlorine molecule. 





254 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


TABLE I 
Vacuum wavenumbers and rotational assignments for the 2780 A band of CCI 








7”— "I 1/2 7— "13/2 

Jo Ra Ri P; Re Riz Q2 P2 Pre 
1.5 35998 . 93 35876.64 35865.60 = 96 

2.5 36000. 10 878.66 866 .57 873.34 35869.63 

3.5 001.29 880.60 867.69 873.80 868.75 

4.5 002.43 35988.98 882.45 868.75 874.38 867.41 

5.5 003.56 987.46 884.36 869.63 874.90 866.57 

6.5 004.71 985.92 886.38 870.63 875.36 865.71 35851.97 
7.5 005.59 984.31 888.43 872.04 875.95 864.88 850.25 
8.5 007.01 982.67 890.49 876.64 864.17 848.49 
9.5 008.12 981.04 892.79 877 .36 863.54 846.81 
10.5 009.30 979.45 894.97 878.13 862.73 844.94 
11.5 010.50 977.79 897.20 878.94 843 .52 
12.5 011.68 976.285 899.41 879.75 841.44 
13.5 012.80 974.75 901.70 880 .60 839.76 
14.5 014.31 973.09 904.02 881.52 837 .97 
15.5 015.48 971.83 906.35 882.45 836.22 
16.5 36044.36 016.69 970.09 908.83 883 .37 834.50 
17.5 046.65 017.99 968.63 911.10 884.36 832.81 
18.5 049.52 019.47 967.17 913.49 885.56 831.09 
19.5 052.24 021.10 965.72 915.92 886.68 829.40 
20.5 054.92 022.23 964.34 918.38 887 .40 827.73 
21.5 057.67 023.65 962.98 920.84 888 .43 826.09 
22.5 060.48 025.10 961.62 923.30 889.54 824.45 
23.5 063.29 26.63 960.33 925.82 890.59 822.83 
24.5 066.09 028.10 959.06 928.35 891.75 821.23 
25.5 069.01 029.72 957.63 930.89 892.76 819.61 
26.5 071.92 031.27 956.59 933.57 894.07 818.05 
27.5 074.87 032.68 936.05 895.23 816.52 
28.5 077.82 034.28 938.59 896.39 814.58 
29.5 080.47 035.91 941.24 897 .64 813.73 
30.5 083.83 037.54 943 .88 898 . 86 811.98 
31.5 086.92 039.16 946.56 900.29 810.51 
32.5 089.98 040.86 949.19 901.41 809.04 
33.5 093.10 042.62 951.89 902.71 807 .64 
34.5 096.26 954.40 904.02 806.20 
35.5 099.49 905 .52 804.82 
36.5 102.63 906 . 66 803.40 
37.5 802.08 
38.5 800.75 
39.5 799.44 
40.5 798.08 
41.5 796 .82 
42.5 795.54 





ROTATIONAL ANALYSIS OF THE 2780 A BAND 


A microdensitometer trace of this band under high resolution is shown in 
Fig. 2. It consists of two subbands separated by approximately 130 cm~', 
with two heads, labelled Q; and Ps degraded to shorter wavelengths. Four 
intense branches can be picked out in each subband, labelled Pi, Q1, Ri, Rei; 
Pio, Po, Qe, Re. The large separation between the subbands, and also the 
presence within each of a fourth chief branch for which necessarily AK # AJ 
implied that the transition was between a Hund’s case (a) and a case (6) 
doublet state. Lines of low J numbering were found to deviate from a para- 
bolic formula,* which implied that the latter state was not purely case (bd), 


*“T ow J deviations” are those caused by the terms in & in equation (3) below. 
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with zero doublet separation, and could not be a ?2 state. The effect of these 
deviations can be most clearly seen in Fig. 2 in the R,; and Qe branches; the 
Q2 branch spacings at low J values are much smaller than those in the R; 
branch, whereas the spacings in these branches should be the same for tran- 
sition from a pure case (6) state. 

Analysis showed that two intense branches were of Q type; hence AA = +1 
for the transition. The electron configuration of CCl demands a *II, ground 
state; electronically similar molecules such as SiF (Eyster 1937), CF (Andrews 
and Barrow 1951), and NO (Schmid 1930) all have ?II,(a) lower states for 
transitions observed under similar conditions. Thus the transition will be 
considered as ?A(b) — 7II,(a), the details of the analysis confirming this 
assignment. All attempts to analyze the band on the basis of a *2 — *II(a) 
transition failed. 


1. Term Values and Expected Spectrum 

Standard analytical formulae are given in the literature (Herzberg 1950 
and references quoted therein) and only certain aspects of these need further 
discussion here. The expected total of 24 branches in the band is reduced 
to 12 if A-type doubling is not separately resolved. All of the main (AK = AJ) 
branches P;, Q1, Ri; P2, Q2, Re were observed as well as the intense (AK # AJ) 
branches R2; and Pj. The remaining weaker branches P21, Qo1; Q12, Riz closely 
paralleled other branches and were not separately resolved, except for a 
number of lines in the R,2 branch that were revealed owing to the low J 
deviations of the Qe branch, as can be seen in Fig. 2. Very few of the P: lines, 
and none of the Q; lines (both head-forming) could be resolved. Thus the 
analysis relied mainly on three branches in each subband (P,, Ri, Ra; Pi, 
Qe, Ro). 

Rotational term values for doublet states intermediate between Hund’s 
cases (a) and (8), relative to a vibronic origin, are given by 


(1) FJ) = Bol(J+3)?F{ (J+3)?+C}}] 


upper and lower signs referring to 7 = 1, (J = K+) andi = 2, (J = K-—3) 
respectively, C? = Y(Y—4)A?/4, with Y = A/B,, A being the spin-orbit 
coupling constant. Additional terms due to A-type doubling and centrifugal 
stretching need not be added to equation (1); from the relation D, ~ 4B,°/w?, 
with approximate values B, = 0.6 cm and w = 860 cm, D,J*? ~ 0.05 cm™! 
even at J = 30 (cf. resolution of ~0.15 cm —'). Also, least squares analysis 
of combination differences with B, as a linear and D, as a higher power co- 
efficient gave values for D, of the order of magnitude of the deviation in 
B,. Magnetic coupling terms y were also ignored for the same reason: inany 
case, they are expected to be negligibly small even in the A state over the 
range of K values used in the analysis. 

For the *II,(a) lower state, with large Y, A, and C, expansion of the square 
root term in (1) as C{1+(J+4)2/C?}? gives 


(2) Fi (J) = By’ /4F By (C’+1/8C") + (Bear) iJ (J +1) 
where the effective rotational constant (Bi), = Bi’ (1F1/2C’]. Quadratic 
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and higher terms in the expansion were ignored since they contribute <0.01 
cm~! for J < 32. Terms independent of J in (2) are not considered in the 
rotational analysis, being included in the vibronic energy. 

For the ?A(6) upper state, with small Y, 4, and C, expansion of the square 
root term in (1) as (J+4){1+(C?/(J+3)*}? gives 


(3) Fi(J) = B Bi J-DU+H)- Ga | 


Fi(J) = Bl oays3/+(z84 te | 


where a = (C’)?/2. Although the original expression (1) is used below in 
determining C’ accurately, (8) was useful in detecting errors in the com- 
bination differences (4) described in the next section. 


2. Lower State Constants 

The analysis began by considering the following two combination differences, 
each of which involved only lower state levels and thus cancelled exactly any 
low J deviations in the upper state. 


R2(J—1)—Q2(J) = 2(Bers)2J 
Ri(J-1)—Pi(J+1) = 4(Bett)s (J +4) 


The assignment was helped by the fact that for incorrect relative numberings 
of branches, terms in a in (3) caused large residual errors in fitting (4) to 
linear functions of J by the method of least squares. The unique numbering 
found for each pair of branches gave the values of (Bii):, (Bus, and 
By = {(By)1+ (BY) 2) /2 shown in Table II. The approximate value obtained 


(4) 


TABLE II 


Rotational constants for the 2780 A band. Errors are 
quoted as standard deviations 


Upper ?A(0) state: By = 0.70533 +0.00062 cm 
ro = 1.63546+0.00073 A 
Doublet splitting: Y’ = —5.87 }o 5 
A’ = —4.14cem71f% | +6.96 cm™ 


Lower 2II,(a) state: (Bit): = 0.68892+0.00101 cm=! 


cit}2 = 0.69506+0.00096 cm= 
By = 0.69199+0.00069 cm™ 
ro = 1.65115+0.00081 cm™ 
Doublet splitting: Y” = 195.03+0.19 
A” = 134.96+0.13 cm™ 
Subband origins: vo) = 36,000.92+0.25 cm™ 


vor, = 35,867,350. 12 cm? 


*Positive values are for a regular, and negative values for an inverted 2A state. 


for Y” from the difference between (Bj): and (Byj)2 was 224+51 cm™, 
consistent with the more precise value of 195.03 cm~ obtained below. 
The doublet splitting A’’(J) is given by the following two relations: 


(5) ANS) = Fi(J)—Fi'(J) = Pi(J)—Pi(J) = Ra(J)—Re(J). 
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The numbering in the P; and R2 branches is already known; it was established 
for the Piz and Re; branches as that for which equations (5) gave identical 
A’ (J) values at each J. From equation (2), 


(6) A"(J) = 2By(C"+1/8C") + { (Bett)2— (Bets) 1} (J+J). 

A least-squares treatment gave 

(7) A” (J) = (133.57+0.013) +J (0.0158 +0.0142) +J? (0.00635 +0.00035). 
The quadratic coefficient agrees with the value (Bi)2—(Bi,)1 = 0.00614 

+0.00139 cm~ obtained from equation (4). The constant term was used to 

obtain the values of Y” and A” in Table II. 


3. Upper State Constants 
The following combination difference gave Bj independently of deviations 


from case (6) in the upper state. 
(8) [Ri(J)—PilJ)]+[{Re(J) — Q2(J)} + {Ra(J—-1) — Q2(J—-1)}] 
=4B;(2J+1). 


The doublet splitting A’ was obtained experimentally using the relation 


(9) A’(N) = [Fi(N)—- Fi(N)] = A”(J)—[Ri(J) - Q2(J)] 

where NV is the upper state value of K; values of A’’(J) were taken from the 
least-squares equation (7). Attempts to fit the result to formulae derived 
from the approximate relation (3) failed, but good results were obtained by 
using the general equation (1), which gives, with NV as the variable replacing J; 


(10) A'(N) = BY(N?+C°)?+{ (W+1)?+C7}?— (2N+1)]. 


A computer program gave C? = 57.98 as showing a minimum residual devia- 
tion of 0.138 cm~! between the calculated and observed values of A’ over the 
range V = 2-34. Values for Y’ and 4’ derived from this are given in Table II. 


4. Determination of Subband Origins 
The origin vo2) of the subband at lower frequency was determined by a 


least-squares fit to the combination relationship; 


(11) Re(J—1) + Pie(J+1) = 2v02)+Bi/2—2(Brir)e 


+2JI(J+1){Bi— (Beis) 2}. 


Since lines close to the other origin voi1) were not resolved, an equation 
similar to (11) would introduce errors of extrapolation. Greater precision 
could be obtained by taking the constant term in (7) as vo(1) —voc2); Compare 
equations (2) and (6). Values obtained are shown in Table II. 


5. Doubling of Lines 

Although A-type doubling was not resolved, detailed examination of lines 
indicated that those in the subband at higher frequency were slightly broader, 
as would be expected for a "II, lower state; doubling in the ?A state would 
be negligibly small. Also, lines in the P12 and Ra; branches, which stand out 
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clearly on opposite edges of the band, showed slight indications of splitting 
into two components. The effect was too small to be resolved properly, and 
cannot be seen in Fig. 2. It did not seem due to A-type doubling, since its 
magnitude was independent of J. Satellite bands do not parallel these branches, 
but the isotopic CCI*’ species could cause additional lines to appear weakly. 
A Zeeman effect caused by magnetic fields in the discharge tube could cause 
doubling of rotational lines, as could quadrupole interactions of the chlorine 
nucleus in one electronic state, but the experimental evidence is too slight to 
decide on this. 


VIBRATIONAL ANALYSIS OF THE BAND SYSTEM 

The rotational analysis shows that the two intense heads of the 2780 A 
band are in the Q; and P»2 branches, with a separation of 138.1 cm~!. Two 
weaker intensity maxima with a separation of 138.5 cm! are ascribed to an 
underlying sequence band labelled (1-1) in Fig. 2. A rotational analysis of the 
more complex band at 2850 A was not performed; nevertheless, similar heads 
can be identified in its high resolution spectrum, and these are shown in 
Fig. 3. Heads can also be identified, though somewhat more tentatively, in 
the two other weak bands of Fig. 1. The frequencies and separations of 
identified Q; and P» heads are listed in Table III. The head separation in a 
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Fic. 3. Microdensitometer profiles of Q, and P2 peaks in the CCI band at 2850 A, under 
high resolution. Peaks marked ‘‘is.”’ are due to CCI’. 
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TABLE III 
Frequencies (cm~') and assignments of band heads and bands (7 = isotopic band) 


Head frequency 
Doublet separation Band frequency 


Sequence Q1 Pe =(0,—P, (Qit+P2)/2 
Av = +1 36836 36700 136 36768 (1-0) 
(2714 A) 
Av =0 35992 .2 35854.1 138.1 35923 .1 (0-0) 
(2780 A) 976.3 837.8 138.5 907.1 (1-1) 
35134.2 35001 .8 132.4 35068 .0 (0-1)z 
Av = -1 130.9 000.1 130.8 65.4 (1-2)z 
(2850 A) 128.2 34995.7 132.5 61.9 (0-1) 
124.6 993.9 130.7 59.2 (1-2) 
Av = -—2 34284 34147 137 34215 (0-2) 
(2920 A) 34246 34098 148 34172 (2-4), etc. 


band will be very similar to the separations of the subband origins, and for 
this vibrational analysis, it is sufficiently accurate to call the mean of the 
head frequencies the band frequency, as shown in the last column. 

Since the internuclear distances in the upper and lower electronic states 
differ by only 1%, we assumed first of all, along with previous workers, that 
the most intense band, that at 2780 A, contained the (0-0) transition and 
its sequences, and the 2850A band contained the (0-1) transition and 
Av = v'—v’’ = —1 sequences. This gives w’’ ~ 860 cm—!. w’ must be of similar 
magnitude if the 2714 A band contains the (1-0) sequence. 

Confirmation comes from the isotope effect. Bands due to the isotopic 
CCI*7 species should appear with appreciable intensity and with nearly the 
same Q,— Pz head separation as the corresponding normal band. The isotopic 
(0-0) band is at almost the same frequency as the normal (0-0) band, and 
buried beneath it, but isotopic (0-1) and (1-2) bands are expected to be 
shifted by +6 cm! from the corresponding normal bands for w’’ ~ 860 cm}. 
This is indeed the shift found in two weaker sets of P2 and Q; band heads, as 
shown in Fig. 3 and Table III. 

If the assignments made here for either the Av = 0 or Av = —1 sequence 
bands are interchanged so as to make the stronger band at higher frequency 
originate at the v’ = 1 rather than the v’ = 0 level, very improbable anhar- 
monicities are found for the lower state vibration. The band at 2714 A can 
now be identified as the first member of the Av = +1 sequence, for the 
vibrational combination relationship predicts 


(1-0) = (1-1) + (0-0) — (0-1) = 36768.4 cm—', 
which agrees with the observed frequency of 36768 cm~! for the band. 
The predicted frequency of the (0-2) transition is 
(0-2) = (1-2)+(0-1)—(1-1) = 34214.0 cm™. 
The 2920 A band is weak and complex, and its features cannot be identified 


with certainty; nevertheless, a band of suitable frequency does occur at 
34215 cm7!. 
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Values so far obtained for the v’’ = 0, 1 levels give the ground state constants 
in Table IV. The first vibrational quantum in the upper state is also found 
to be 845 cm“. 


TABLE IV 


Vibrational constants 


2A 211 
We 860 875.1 cm™ 
Were ~7.5 7.0 cm 
D, ~24,700 27,300 cm=!* 


*Calculated from Dg = we?/4weXe. 


The remaining peaks in the 2714 A band at 34246 and 34098 cm! have 
roughly the correct separation for a Q; and P2 head. They do not belong to 
the (1-3) band, since they are more intense than the heads of the (0-2) 
band, and also this assignment would give 887 cm~! for the third vibrational 
quantum in the lower state, implying negative anharmonicity. However, if 
wx, is given the plausible value of ~7} cm™', it can be shown that heads 
in successive bands in the Av = —2 sequence will themselves converge at 
the observed head frequencies, and that this type of convergence will not 
occur for successive bands in other observed sequences. This accounts for 
the relatively high intensity of the two peaks, and is consistent with the rest 
of the vibrational scheme. 


DISCUSSION 


In their rotational analysis, Kuzyakov and Tatevsky (19596) ignored low 
J deviations of lines from parabolic formulae. This led to a wrong numbering 
of lines in the R2 and Q2 branches upon which their analysis is based. Also, these 
deviations invalidate the band head-origin formulae used in their vibrational 
analysis (1959a). 

There remains the possibility that the transition is *II(6) — ?II(a), with 
A-type doubling of the branches remaining too small for observation. Addi- 
tional rotational lines expected near band origins for such a transition were 
not identifiable, although the resolution was not sufficient for this to be 
conclusive. The chief argument in favor of a *A upper state is that the Q 
branches in the bands have high intensity. In Fig. 2, the Q2 branch, which is 
distinct from its satellite Ri. branch at low J, shows strongly between the 
wider-spaced lines of the R: branch. The high intensity of the head in the 
Q, branch cannot be ascribed to the accompanying satellite P2; branch, since 
it can be calculated that the latter would form a head 1 cm! to lower fre- 
quencies, and such a head does not appear with appreciable intensity. 

The upper state is identified as 7A; rather than ?A,. The lines Q2(1.5) and 
P,(2.5), both originating from the J’ = 1.5 upper state level, can be assigned 
in Fig. 2 in the correct positions for an inverted upper state. The alternative 
values of A’ and Y’ for a regular state given in Table II would cause both 
these lines to appear 11.1 cm~! to lower frequencies, where they are not 
observed. 
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The lower state of the transition is almost certain to be the ground state 
of the radical, arising from the . . . r4e?x electronic configuration. The ground 
state bond length of 1.65 A lies between the approximate values of 1.72 A 
and 1.63 A for C—Cl bonds adjacent to double and triple bonds respectively 
(Bowen 1958). Similarly, Andrews and Barrow (1951) found that the ground 
state internuclear distance for the CF radical was 1.27 A, slightly less than 
for double-bonded CF groups (~1.32 A). 

The ground state dissociation energy (Do = w?/4woxo) of 3.34 ev is less 
than the CF value of 4.96ev as would be expected. The excited state 
Dy ~ 3 ev gives approximately 4 ev as the difference in energies between the 
dissociation products in the two states, but this value results from rather 
drastic extrapolations. Assuming, as Andrews and Barrow have for CF, that 
the lower state gives normal C(*P) and halogen (?P) atoms, then the upper 
2A state found here could correlate with either; firstly, the same atomic states 
C(8P)+Cl(?P) as in the ground state; secondly, the states C(!D2)+Cl(?P); 
or thirdly, highly excited states of both the C and Cl atoms, with the possi- 
bilities here too numerous to specify. It is improbable that the extrapolated 
dissociation energies are sufficiently in error to make their difference equal to 
zero, but the difference of 1.2 ev required by the second alternative makes 
this a more likely possibility. 

Finally, although the shorter bond length in the excited state is associated 
with a lower vibrational frequency than in the ground state, the differences 
in frequency and bond length are only 1-2%. Badger’s rule relating bond 
length and frequency cannot be expected to apply when such small differences 
are involved. Interestingly, Eyster (1937) specifically mentions the poor 
results obtained when applying Badger’s rule to doublet states of the SiF 


molecule. 
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ON THE THERMAL EXPANSION OF METALS AT 
LOW TEMPERATURES! 


G. K. Horton 


ABSTRACT 


A theory is developed which correlates the thermal expansion of crystals to 
the anharmonicity introduced into Born’s lattice dynamics by allowing the force 
constants of the crystal to vary with volume. This is achieved by identifying the 
force constants with the elastic constants of the crystal by the method of long 
waves. It is then assumed that it is primarily the volume dependence of the 
elastic constants that give rise to their temperature variation. A central force 
nearest and next-nearest neighbor force model analogous to Leighton’s is applied 
to copper. The values of the lattice thermal expansion coefficient and of Griinei- 
sen’s parameter are given as a function of the temperature and found to agree 
quite well with the latest experimental results. It is pointed out that the descrip- 
tion of the interionic potential in metals by a two-body central force is certainly 
a serious oversimplification and that the theory is likely to be more realistic for, 
say, the ideal inert solid gases, as soon as the experimental data becomes avail- 


able. 


I. INTRODUCTION 


The lattice dynamics of Born (Born and Huang 1954) has met with con- 
siderable success in those effects that depend primarily on the harmonic 
contribution to the free energy. Among these I need only refer to the theory 
of specific heats where Born’s theory is able to account for the considerable 
body of experimental information with reasonable success (Blackman 1955). 

Of the various phenomena that depend on the anharmonic terms in the 
free energy, I wish to discuss the thermal expansion of crystals. This seems 
appropriate in view of the recent experimental work in this field. In an ideally 
harmonic crystal, e.g. one in which all pairs of atoms attract or repel each 
other with forces proportional to the distance between them, there would be 
no thermal expansion. No complete theory of lattice vibrations that takes 
proper account of anharmonic effects has vet been carried through satis- 
factorily for a three-dimensional dielectric crystal, much less for a three- 
dimensional metal. I would like to emphasize that this paper makes no con- 
tribution to this problem. 

As has been discussed by Salter (19555), Born’s lattice dynamics (Begbie 
and Born 1947) is only quasi-harmonic and, in a certain sense, has anharmonic 
effects built into it. In this paper I present an investigation of the extent 
to which a suitably formulated quasi-harmonic theory can account for thermal 
expansion of crystals. In so doing, I would like to stress the exploratory nature 
of the discussion. Some version of Born’s theory, in which the interionic forces 
are represented by one general potential function, with due reference to the 
influence of the conduction electrons (Bhatia 1955; de Launay 1953) on this 
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potential, may, in principle, be able to account for the low temperature heat 
capacity of crystals. This may not be true of the thermal expansion and so 
we probably cannot expect quantitative agreement with experiment. 

The theory of the equation of state of crystals is largely due to Mie (1903) 
and Griineisen (1912). This theory, in its usual form (Griineisen 1926), is 
based on the well-known Debye approximation, leading to Griineisen’s rule 
that the thermal expansion coefficient and the specific heat of a crystal are 
proportional. Recently Barron (1955a) has studied the effect of replacing the 
continuum model of Debye by the dynamical version of the nearest neighbor 
central force theory of Born and von Karman (1912). Barron found that in 
his theory, Griineisen’s parameter y was no longer a constant at intermediate 
temperatures, though it did approach (different) constant values at high 
(y.) and low (yo) temperatures. The high temperature value y,, depended 
on the logarithmic volume derivative of the one force constant that enters into 
the lattice dynamical model used by him. It is left as an arbitrary parameter 
(which may still vary with temperature through its volume dependence) that 
has to be found for a given crystal by identifying the experimental and 
theoretical y,,. The temperature dependence of y in the intermediate region 
was then a function that depended only slightly on the properties of particular 
crystals. Using a Mie—Lennard-—Jones potential, he also estimated the influence 
of more distant neighbors: this influence turned out to be quite marked. 

In this paper I would like to report on an attempt to use the nearest and 
next-nearest neighbor central force theory (Leighton 1948) in its dynamical 
form together with the experimentally measured elastic constants. It is in this 
way that we introduce the major anharmonic effects, including that of the 
zero point energy, into the theory. The theory presented in this paper contains 
no adjustable parameters and yields results that can be compared directly 
with the experimental results for the crystal whose elastic constants are being 
used. It is applied to a face-centered cubic metal Cu, because, when the 
numerical work was started, it was only for this crystal that reliable values 
for the expansion coefficient, the heat capacity, and the elastic constants 
were available. I would like to emphasize that this does mot mean that I 
believe that the interionic forces in a metal can be described by a central, 
two-body short range force (cf. Peierls 1955). An application to the heavier 
inert gas solids would clearly be more realistic when the experimental results 
are available with sufficient accuracy. The choice of the nearest and next- 
nearest neighbor central force theory was indicated by two considerations. 
The first is the reasonable success that this theory enjoys in accounting for 
the specific heat of Cu and, to a lesser degree, of other noble metals (Horton 
and Schiff 1959). The second arises from the fact that this theory is the 
simplest in common use that still permits direct comparison with experiment. 
The assumption made in the theory that each “‘shell’’ of neighbors is separately 
in equilibrium is made for simplicity in this exploratory study. In future work, 
especially when the zero point energy plays an important role, it should not 
be made (Barron 19556). 
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II. THEORY 


The Mie-Griineisen equation of state is given by 
; 3N 
(1) PV = DO+kT DI ysP(x}), 
j=l 
where ® is the potential energy of the crystal with the particles at rest in 
their equilibrium positions, x, = hw,/kT, P(x) = x/(e7—1), 
(2) Vi = —Qln Wi, 


QZ = V(0/dV), and all other symbols have their usual meaning. w, and ® 
are functions of the volume only. From the usual thermodynamic relations 
we find from (1) that 


(2) _ Cs 
(3) eae 





Here a = (1/a)[(0a/9T)p] is the linear thermal expansion coefficient, and 
x = —(1/V)[(@V/aP)7] is the isothermal compressibility, a is the nearest 
neighbor distance and 


3N 
(4) = 2X VsE(x5)/Co, 


with E(x) = x%e7/(e7—1)?. 

The dynamical matrix for our central force theory may be found in Salter 
(19552). The summations occurring in (4) were replaced by integrals in the 
usual manner. These integrals were then evaluated by using the modified 
Houston spherical five-term integration method developed by Horton and 
Schiff (1959). The error introduced by the use of this method does not amount 
to more than about one per cent over the whole temperature range considered. 
To use the integration formula we must solve the secular equation in reci- 
procal lattice space and insert the resultant dispersion law into the integrals 
obtained from (4). I will carry through the procedure in the A, or (100), 
direction as an example of the general method. One finds that 


= [2(261+61) + (82+52)p](2—p), 


(5) (5) a long 


ll 


(6) (a) cues + [(48:+81)+B2p](2—p), (twice) ; 


p = (1+ cos ka/+/2), p is the density, and a is the nearest neighbor distance. 
_ pa*[ 1 a : 
(7) Bi = on } OP Jone’ 


3 
_ sat) ; 
(8) n- 2/4 ae 
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corresponding expressions for 82 and 62 can be found by replacing a by +/2a. 
In the static theory we put the #’s (proportional to the first derivative of the 
interionic potential) equal to zero and, using the usual Fuchs—Fine method, 


find 


61 = Ca, 


(9) 
be 


(C11 — C12 — Ca) /2. 


Because of the form of (2), we cannot use this procedure here. Clearly, to 
find the y;, all volume differentiations must be carried out before the §’s are 
equated to zero. From (2) we note that 


(10) 1s = — G2 Pao 
and that 9(+) = a . 
a°p p 
d 
i 9,=a—. 
with - a 


From (7) and (8), and noting the volume dependence of the density, it is 
easy to show that 


(11) DB; = (6:—8:). (¢ = 1, 2) 


Next we must consider 2,5; (i = 1, 2). This expression can be evaluated 
directly if one assumes a definite function for W. Our problem is to evaluate 
it subject to (9). We do this with the assumption that 


Ne) 
dT/\da/’ 
1 dé, 


(12) D8; 


the pressure being constant. 

This assumption may be regarded as equivalent to the statement that it 
is assumed that if the elastic constants were measured at constant volume as 
the temperature varies, they would essentially be temperature independent. 
We assume that this is a reasonable approximation well below the characteristic 
temperature. Finally, following Barron (1955a), we will ignore the volume 
dependence of ~ in (5) and (6). After the appropriate differentiations corre- 
sponding to (10) have been carried out, the 8;’s can then be equated to 
zero; they are likely to be very small. In this way we obtain the +, corre- 
sponding to (5) and (6). They are given by 


(13) (v3) long ~ 1 25,+5op 7 26: +50p ’ 


(14) (y;)* trans ~ a2 i ’ (twice). 


6; a 6, 


| 


GF MOIR NRAIE I so 


oe 


i 





are 





HORTON: THERMAL EXPANSION OF METALS " 267 


The dot over the 4’s indicates a derivative with respect to the temperature. 
The results for the other four directions (110), (111), (210), and (211), may 
be derived in a similar fashion. The values of the elastic constants and their 
temperature derivatives are taken directly from the experimental results. For 
the temperature range 0° K to half the Debye characteristic temperature it 
is convenient to fit the published tables by a quadratic polynomial in T by 
the usual methods. At very low temperatures this procedure breaks down. The 
breakdown was traced to the fact that for T < 6/10 it becomes important 
to represent the elastic constants by the correct theoretical expression (Born 
and Huang 1954). In a second calculation I, therefore, wrote 


(15) Ciyy = a—bT'+cT® (7j = 11, 12, 44). 


The results of the two calculations can then be used in the respective tem- 
perature ranges in which they are valid. 

At high temperatures finally we can use the moments method in the form 
first given by Barron (1955a) and Horton and Schiff (1959). We write 


3 


(16) E(x) = YS Anx™, 


n=0 
with 
(17) Ao=1, Ae = —0.08320, A, = 0.003977, As = —0.000106. 
Equation (17) has been obtained from (16) by a fitting procedure accurate 


to about 0.1% in the range 0 < x < 3. All the higher 4’s are zero. On inserting 
(16) into (4) we find 


3 3 
(18) > » A on'V2nl2n fe A nC ons 


n= 


where 
(19) Con = Mon(h/kT)™, 
where 
1 2 
(20) Hn = ay 2, @j'5 
and 
(21) Yn = 5 FIN (n > 0). 


The case n = 0, which gives y,,, requires special treatment. Clearly 


1 3N 
(22) Yo = ay dy Ye 


j=l 


Equation (22) may be evaluated directly by converting it into an integral. 
Indeed the integral is just the numerator of (4) with the Einstein function 
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replaced by unity. Equations (13) and (14) are samples of the integrand that 
can again be evaluated by the modified Houston method. 

For higher terms in (18) we note that after carrying out the indicated 
differentiations and with 


. ae 
(23) 6.. = k (/ 3 Ha 


(24) n= a4 (Ait in/2)/(4-+54/2) | : 


(25) a 5 (0,/T)* 


Higher terms in the expansion are similar and will not be reproduced here 
explicitly. 

Both in the high temperature limit and at lower temperatures we see that 
the results of our analysis of (3) and (4) can therefore be put in the form 


(26) . a@=x(q+r/a), 
where 
(27) x = x/8V. 


The g and r are quantities that come from the numerator of (4) corre- 
sponding respectively to the portion of y, in (13), (14), and (24) that is 
independent of a and that multiplies 1/a. Equation (26) is quadratic in a 
with solution 


(28) a = x{g/2+V7[(q/2)?+r/x}}. 


It can be shown that we must choose the positive sign of the square root 
in (28) as it corresponds to a higher entropy state, i.e., one with a lower free 
energy (Varley 1956). We obtain Griineisen’s parameter as a function of 
temperature from (3), (26), and (28). The result is 


(29) y = a/xC,. 


All the numerical integration was done on the University of Toronto 
Ferranti I digital computer and on an LGP 30 Royal McBee computer. It is 
a pleasure to thank Canadian National Telegraphs for making a direct teletype 
line available to the University of Toronto Computation Center and to thank 
the National Research Council for free time on Ferut. The integrations were 
performed at intervals of 10° K by using the wave number integration pro- 
cedure, combined with Houston’s spherical five-term formula and Simpson’s 
rule with 24 intervals at 10° K and 12 intervals at higher temperatures. The 
detailed procedure is discussed in Horton and Schiff (1959). 


III. DISCUSSION OF THE RESULTS FOR COPPER 


The elastic constants for Cu were taken from Gaffney and Overton (1955). 
I used a = 2.551X10-8cm and p = 8.995 g/cm*. The experimental results 
are those of Rubin, Altman, and Johnston (1954), Simmons and Balluffi 
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(1957), and White (1960). As the results of the capacitor dilatometer method 
as used by Bijl and Pullan (1955) seem to contain certain systematic errors 
at low temperatures, I confine the discussion to the results of the former 
authors. The thermal expansion coefficient as a function of temperature is 
given in Fig. 1 while the temperature dependence of Griineisen’s parameter 
is shown in Fig. 2. 


T 
x 
° 
° 
2 
= THEORY (WITH QUADRATIC 
x POLYNOMIAL) 


THEORY 
(WITH oi 





Oo 20 40 60 80 100 120 
——~s> T °K 


Fic. 1. The linear thermal expansion coefficient of Cu as a function of the temperature. The 
theory (with quadratic polynomial) refers to the description of the temperature dependence 
of the elastic constants in terms of a quadratic polynomial in 7. Equation (15) incorporates 
the correct 7‘ low temperature dependence of the c;;. This curve is only valid below to 30° K. 
The experimental curve was obtained from the measurements of Rubin, Altman, and Johnston 
(1954), Simmons and Balluffi (1957), and White (1960). 


It is clear that the theoretical results at higher temperatures lie about 8% 
above the experimental results. At lower temperatures the values of y become 
rather uncertain theoretically as I am trying to calculate the ratio of two 
quantities that separately tend to zero as 7%. The increase in y below 40° K 
is clearly due to this cause and hence probably spurious. In this region the 
experimental results also become increasingly less accurate. It is neverthe- 
less clear that the dashed curves obtained by using (15) to represent the 
temperature dependence of the elastic constants leads, as might be expected 
at very low temperatures, to more reasonable results. At 30° K neither of 
the two theoretical curves are reliable and an interpolation between the two 
must be made. In view of the crude theoretical model I feel the agreement 
between theory and experiment is surprisingly good. 

In conclusion I would like to enumerate briefly what seem to me to be 
the most serious defects of the theory presented above. With regard to the 
theory in general, it is firstly phenomenological in the sense that I have 
simply related two phenomena that depend essentially on the anharmonicity 
of the lattice vibrations, namely the temperature variation of the elastic 
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THEORY (WITH 15) 
THEORY (WITH QUADRATIC POLYNOMIAL) 
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EXPERIMENT 





Oo 20 40 60 80 100 120 140 160 180 


Fic. 2. The temperature dependence of Griineisen’s parameter y. The theory (with quad- 
ratic polynomial) refers to the description of the temperature dependence of the elastic con- 
stants in terms of a quadratic polynomial in 7. Equation (15) incorporates the correct T* 
low temperature dependence of the c;;. This curve is only valid up to 30° K. The experimental 
curve was obtained from the measurements of Rubin, Altman, and Johnson (1954), Simmons 
and Balluffi (1957), and White (1960). 


constants and the thermal expansion. A proper theory (Ludwig 1958) would 
contain an explanation of both phenomena. Secondly the temperature varia- 
tion of the elastic constants may not be due to the expansion of the solid 
only. It is likely that any additional corrections are small, however, because 
the elastic constants of a number of crystals, at constant volume, appear to be 
reasonably constant with temperature. Thirdly the use of the modified Houston 
method may have introduced errors of one or two per cent especially at low 
temperatures. Fourthly the assumption of nearest and next-nearest neighbor 
central forces is a considerable oversimplification. There is no difficulty, in 
principle, in removing these last assumptions. As far as the application to 
copper is concerned, this was done purely for the purpose of illustrating the 
kind of behavior of a and y that is predicted. Throughout we have ignored 
the explicit contribution of the conduction electrons to the thermal expansion 
which is certainly a good approximation for Cu but may be questionable in, 
say, Al (Varley 1956). 
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REFRACTION AND DIFFRACTION OF PULSES! 
W. E. WILLIAMS 


ABSTRACT 


It is shown that the solution for the reflection and refraction of a plane pulse at 
a plane interface may be obtained, by Fourier synthesis, from the correspondin 
solution for a time-harmonic plane wave even when the angle of incidence is suc 
that total reflection occurs. The application of the Fourier superposition method 
to problems involving both refraction and diffraction is also discussed. 


1. INTRODUCTION 

It is a standard method in the theory of scattering of waves by obstacles to 
obtain the solution for incident plane pulses of arbitrary shape by Fourier 
superposition of the solution for an incident time-harmonic plane wave (for 
example, Wait 1957). It has, however, been pointed out (Friedlander 1948) 
that, if the angle of incidence is such that total reflection occurs, this approach 
breaks down even in the simplest problem of reflection and refraction at an 
infinite plane interface. The object of the present work is to indicate the method 
by which the Fourier synthesis approach may be employed to solve this type 
of problem for all values of the angle of incidence. If a time variation 
exp—iwt is assumed then it is shown that the reflection and transmission 
coefficients for the case of total reflection depend on the sign of w. The results 
for arbitrary plane pulses can then be deduced from known results in the theory 
of Fourier integrals. The results are in agreement with these obtained by 
Friedlander (1948), who obtained the solution in an entirely different manner. 
The apparent breakdown of the ordinary Fourier superposition method is due 
entirely to the dependence of the reflection and transmission coefficients on 
sign w. This prevents a direct and simple relation being obtained between the 
incident, reflected, and transmitted waves. The dependence on sign w is due 
to the fact that each separate Fourier component, for total reflection, has to 
be evanescent away from the interface and thus the appropriate form of 
each component will depend on the sign of w. 

The problem of diffraction of a plane pulse by a half plane at the interface 
between two media is also discussed. A problem of this type has been considered 
by Papadopoulos (1959) when the incident wave is a step function. Papa- 
dopoulos employs arguments of dynamic similarity in order to effect a solution; 
this type of argument does not seem to be capable of generalization to plane 
pulses of arbitrary shape. 

The problem for the case of an incident time-harmonic plane wave is of the 
type which may be solved by the Wiener—Hopf method. It is shown that the 
results for reflection at a plane interface and the method of Fourier synthesis 
enable the solution for a plane pulse to be deduced directly from the time- 
harmonic solution for all values of the angle of incidence. 


1Manuscript received August 30, 1960. 
Contribution from the Department of Applied Mathematics, The University, Liverpool 3, 


England. 
Can. J. Phys. Vol. 39 (1961) 
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2. REFRACTION OF A PLANE PULSE 
The mathematical problem considered is the solution of 


ee. O46. 1.26 


(1) ax! t ay ae ee y>0, 
(2) ag , HG _ 1 db 4 aa 
ax? iy a: no! ¥ ’ 


subject to the conditions 


bidyao+ = dbadyno-, 


2) aC) 
ai(28 on “ Oy / y—0- R 


where dj, @2, 01, be, C1, C2 are all real constants. It is assumed that there is 
incident in the region y>0 the plane wave 


$ -s(: asin etna) 
= spor Snead aetna 


a 


where ¢ denotes the time. In the present work we are interested in the case 
when the angle of incidence, 0, exceeds the critical angle, i.e. @,;<az sin 6. 
The total solution is given by 


got i, y>0d, 
@ — 
2, x < 0, 


where ¢; and ¢2 are solutions of equations | and 2 respectively and represent, 
in some sense, waves going out from the plane y = 0. 

The special case f(¢) = exp —iwf will be considered first, the solution for 
arbitrary f will then be obtained by Fourier superposition. The appropriate 
forms for ¢; and @2 are 


hep pp ee 


oi = 
ay, 
ane x. 
¢o: = Bexp — tat —sin 6+AYV(, 
a, J 
where 
* = A a s 
ao a, 


and A and B are constants. \ is purely imaginary and thus in order to obtain 
outgoing waves we must have A = (78/a:) sgn w where 6 = +/(sin? 6 — (a?/ 
a;)) and sgn w denotes the sign of w. 

The formal solution of the boundary value problem is elementary and we 
obtain 
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= b3ci cos’ 6@—bic28°+21b1b2¢1628 cos 0 Sgn w, 


Alw Als 


= 2b1¢, cos 0(b2c; cos 0+7b1628 sgn w), 


where 
C = (bici cos” 6+bic28"). 
We are now in a position to solve the problem for an arbitrary plane wave, 
it will be assumed that f(¢) possesses a Fourier transform F(w) in the ordinary 
sense, i.e. 


1 ° aia 
10) = Sep J Pere “ae. 
The reflected wave is thus 
(3)  (b3ci cos” 0—bic38") Cf (¢) +2ib1b2c1¢28 cos 15 - f F(w)sgn we “dw, 


where ¢ = ¢—[(x sin 6 — y cos 6)/a,]. From equations 5.1.8 and 5.1.11 of 
Titchmarsh (1937) it is seen that the above integral is equal to 


bp fli, 


rg ~I-f 


where P denotes the principal value. The transmitted wave is 


(4) bic, cos 0C - F(w) [bec cos 0+7b1628 sgn w] 
Tv “ca 
exp—ta( ~ sin 9+ ty sgn w) ds 
a, a 


From equations 5.13, 5.15, 5.3.5, 5.3.6 of Titchmarsh (1937) we have that 


(5) ‘< F(w)exp—olif+u sgn wldw = z ie : 
(6) I: F(w)exp—ofit+p sgn w)sgn wdw = a 7 fa pioe 


Equations 3, 4, 5, and 6 represent the formal solution of the problem. By a 
suitable choice of the constants 0,, bo, ¢1, C2 the solution for the reflection of 
acoustic, electromagnetic, or elastic plane waves may now be obtained. The 
case b; = b2 has been treated by Friedlander and, for this case, the present 
solution agrees with that of Friedlander, who used a completely different 
approach. The above analysis is strictly formal and is clearly rigorous when 
f and its derivatives possess Fourier transforms in the ordinary sense. A similar 
type of restriction was also imposed by Friedlander. Use of generalized function 
theory (Lighthill 1958) enables the analysis to be applied to other functions, 
for example, the Heaviside unit function. 


3. DIFFRACTION PROBLEMS 


We now consider briefly the application of the results of Section 2 to prob- 
lems involving both diffraction and refraction. Ti. zeneral boundary value 
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problem considered is the determination of a solution ¢ of equations 1 and 2, 
satisfying the boundary conditions of Section 2 on y = 0, x < 0 and with 
¢ or 0¢/dy vanishing on y = 0, x > 0. Boundary value problems of this type 
occur in the diffraction of acoustic waves by either a perfectly soft or per- 
fectly rigid semi-infinite screen placed at the interface between two different 
media. A similar type of boundary value problem occurs in electromagnetic 
diffraction (Clemmow 1953). 

The solution ¢ for an incident plane wave exp —iw{t — [(x sin 86+ cos @)/a,]} 
is easily obtained by Wiener—Hopf methods. ¢ may be written as a sum of the 
plane wave terms of geometrical optics and a diffracted term. We consider 
first the diffracted term and examine the form of this in y > 0, the region 
y < 0 may be treated similarly. From the results of Clemmow (1953) it can 
be seen that the diffracted field has the form 


J ef teria) cosh 9 F(v,)d2, 
0 


where x = rcosy, y = rsiny, and F is independent of w. Fourier super- 
position thus shows immediately that the diffracted field due to the incident 


plane pulse ¢o is 
J hen cosh 0) Flop yao 
0 a) 


The geometrical optics term due to reflection at y = 0, x < 0, is an ordinary 
plane wave of amplitude +1 depending on the boundary conditions imposed. 
The corresponding term for the plane pulse problem is thus the reflected 
pulse term of geometrical optics. It thus only remains to consider the geo- 
metrical optics terms due to reflection at y = 0, x > 0. These have precisely 
the form of ¢; and ¢2 and may be treated in the same manner. Thus, in the 
time-dependent diffraction problem, we obtain the reflected and refracted 
pulses of geometrical optics for a; > a2 sin @ and the solution for a; < a2 sin @ 
is given by equations 3-6. 

The above analysis is quite general and the particular form of F(v,~) will 
depend on the boundary conditions, and obtaining it will involve a complicated 
Wiener—Hopf factorization. The object of the present work, however, is not to 
consider a particular case in detail but merely to indicate how the Fourier 
synthesis technique may be applied, even under conditions of total reflection, 
to obtain the solution of time-dependent problems from the time-harmonic 
solutions. It is of interest to note that, when f(¢) is the unit step function 
H(¢), the diffracted part of 0¢/dt has the particularly simple form H(a,t — r) 
F(cosh~!(a,t/r),w); this is precisely the form of solution obtained from argu- 
ments of dynamic similarity and conical flow methods. 
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EXCITED STATES OF O'8 STUDIED BY THE 
REACTION H*(O}%,py)O'8! 


A. E. LITHERLAND,? R. BATCHELOR,’ A. J. FERGUSON, AND H. E. Gove 


ABSTRACT 


Gamma rays from the excited states of O at 3.63 and 3.92 Mev have been 
observed using the reaction H*(O}8, py)O" at an incident O!* energy.of 14 Mev. 
Both states were observed to emit gamma rays to the 1.98-Mev 2+ first excited 
state of O'8. No evidence for crossover transitions was found and in each case 
the crossover transition was estimated to be <15% of the cascade transition. 
Angular correlations of the gamma rays were obtained and these strongly sup- 
port an assignment of spin 0 to the 3.63-Mev state and a spin of 2 for the 3.92- 
Mev state. These assignments have been confirmed by a recent experiment on 
the O'8(H3,p)O"8 reaction which gives the assignments 0+ and 2+ for these two 
states. Thus the states at 3.55, 3.63, and 3.92 Mev form a triplet with assignments 
4+, 0+, and 2+ which strongly resembles the vibrational spectra found in 
heavier nuclei. However, the measured angular correlations of the gamma rays 
from the 3.92-Mev state show only a small admixture of electric quadrupole 
in the 1.94-Mev gamma ray with relative amplitude +0.1+.1. A lower limit 
of ~10-” seconds on the lifetime of the 3.63-Mev state was obtained from the 
absence of a doppler shift of the 1.65-Mev cascade gamma ray. 


INTRODUCTION 


The excited states of O'8 are of particular interest because of the possibility 
of describing them as the result of the interaction of two neutrons outside the 
closed shell nucleus O". Intermediate coupling calculations by Redlich (1954) 
and by Elliott and Flowers (1955) have been successful in explaining some 
features of the O'8 energy level spectrum. Both calculations predict similar 
sets of low-lying states for the two neutrons in the s,d shell: a 0+ ground 
state and low-lying excited states of 0+, 2+(two), and 4+. Higher excited 
states of spin 0+, 1+(two), 2+(three), 3+(two), and 4+(one) are also 
predicted. 

The wave functions derived by Redlich for O'S have been successfully used 
in a comparison with theory of the ratios of the reduced stripping widths 
determined from the O"(d,p)O'8 reaction by Bilaniuk and Hough (1957). 
Bilaniuk and Hough concluded that the 1.98- and 3.55-Mev states probably 
had spins 2 and 4. The parities of these states were found to be even because 
of the prominent / = 2 stripping patterns observed. Subsequently, in a study 
of the C'(a,y)O!8 reaction by Gove and Litherland (1959), the 1.98- and 
3.55-Mev states were shown to have spins of 2 and 4 as conjectured by Bilaniuk 
and Hough. The spin of the 1.98-Mev state has also been established by the 
work of Phillips (1958) on the C(a,y)O'® reaction. Recently Silbert and 
Jarmie (1959) discovered another level in O'8 at 3.65 Mev between the known 
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levels at 3.55 Mev and 3.93 Mev (Ajzenberg-Selove and Lauritsen 1959). 
This new level was strongly excited by the O1%(t,p)O'8 reaction. The available 
information on the low-lying levels in O'8 is summarized in Table I, which 


TABLE I 


The energies of the low-lying states of O are 
tabulated together with the known spins and 
parities. The energy measurements are (a) from 
the tabulation of Ajzenberg-Selove and Laurit- 
sen (1959) and Silbert and Jarmie (1959) and 
(b) Jaffe et al. (1960). The spins and parity 
assignments are from several sources as discussed 
in the text 


Energy (Mev-kev) 


Spin and parity (a) (b) 
0+ 0 0 
2+ 1.982+4 1.97945 
44 3.55+20 3.55245 
0+ 3.65420 3.63445 
2+ 3.93440 3.915+5 


4.45+50 4.448+5 


includes the most recent values of the energies which were measured by Jaffe 
et al. (1960). The energy values quoted by Jaffe et al. (1960) are used in 
this paper. 

Since this experiment was completed some preliminary results (Jaffe et al. 
1960) have become available on the angular distributions of the proton groups 
from the O'*(¢,p)O'8 reaction. Using 5.5-Mev tritons the angular distributions 
of the protons to the 3.63- and 3.92-Mev states in O'8 suggested the assign- 
ments of 0+ and 2+ for these states. These angular distributions were inter- 
preted on a recent double stripping theory by Newns (1960). 

This paper describes the measurement of gamma-ray angular correlations 
from the 3.63- and 3.92-Mev states and their interpretation in terms of the 
properties of the states. In Appendix I it is shown that the discovery of a 
0+ state 80 kev above the 3.55-Mev state studied by Gove and Litherland 
(1959) does not alter the spin assignments that they made. In Appendix II 
an estimate is given of the lifetime of the 3.63-Mev state based on the absence 
of a doppler shift of the 1.65-Mev reaction gamma ray. It is also pointed out 
that the type of reaction H*(O}$,p)O'’, used in this work is potentially of 
great value in such lifetime measurements. 


EXPERIMENTAL APPARATUS 


To produce O'8 in its various excited states a beam of O" ions was used to 
bombard a thick tritium target absorbed in a zirconium layer deposited on a 
thick copper backing. O'* ions with three, four, and five electrons removed are 
readily produced by the Chalk River tandem accelerator. Beams of these ions 
have previously been used to calibrate the momentum scale of the 90° deflecting 
magnet of the accelerator by means of the reaction H?(Oj$,2)F" (Gove et al. 
1958). Currents of approximately 0.4 microampere of O/{ ions were used for 
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most of the experiments described in this paper. The energy of the O/§ ions 
used was 14 Mev which is approximately the highest energy that can be used 
before gamma-ray cascades in F'8, from the competing H*(Oj§,2)F'** reaction, 
interfere with the measurement of the gamma rays from the 3.63- and 3.92-Mev 
states in O'8, The Q-value of the H*(Oj§,”)F'® reaction is 1.28 Mev and at 
14-Mev bombarding enerzy the energy in the center of mass system is 2.21 
Mev. Consequently only states in F'* up to the 3.35-Mev state can be excited. 

To produce the 14-Mev OS ions, the tandem accelerator accelerates negative 
oxygen ions to the central positive terminal of the machine. Several electrons 
are then removed in a gas stripping canal and the resulting positive ions are 
further accelerated. Charge states other than O/§ have been observed and, 
for example, currents of approximately 0.15 microampere of 36-Mev O}* have 
been used in Coulomb excitation experiments. Small beams of O¢f have also 
been obtained. cs ions are available from the tandem machine when COs is 
used in the ion source. 

After a deflection through 90° by the beam analyzing magnet the O/f ions 
then travelled for about 40 feet into a well-shielded room. Two alternating 
gradient magnetic lenses similar to the type described by Bromley and Bruner 
(1954) were used to focus the beam over this distance. The well-shielded room 
contained an angular correlation table which supported two Nal(Tl) gamma- 
ray detecting crystals: one 5in. diameter by 4 in. long crystal mounted so 
that it could rotate about the target in a horizontal plane, and one 5 in. 
diameter by 6 in. long crystal mounted on another support which permitted 
it to be moved to any point on a hemisphere centered on the target and lying 
on and above the horizontal plane of the other counter. At all such points 
the axis of the gamma-ray counter intersected the target center. The target 
chamber was constructed so that it had nearly equal gamma-ray absorption 
properties for all accessible polar angles 6 and ¢ of the gamma-ray counters. 
The distance from the target to the front faces of the crystals was 6.2 inches. 

The electronic system used to study the pulses from the Dumont 6365 
multipliers which were optically coupled to the NaI(T1) crystals was a con- 
ventional ‘‘fast-slow”’ system (Bell, Graham, and Petch 1952) with a resolving 
time, 27, of 50 nanoseconds. This system permitted a coincident gamma-ray 
pulse spectrum to be displayed on one 100-channel transistorized pulse-height 
analyzer and a direct pulse spectrum to be displayed on another 100-channel 
pulse-height analyzer for monitoring purposes. 

For a measurement on the gamma rays in coincidence with the neutrons 
from the reaction H*(O/$,my)F'8 an NE 212* organic liquid scintillator, 3 in. 
diameter by 3 in. long was used in conjunction with the 5in. diameter by 
4 in. long Nal(T1) crystal. To ensure that coincidences between gamma rays 
were not observed a pulse-shape discrimination circuit of the type described 
by Owen (1959) was used to eliminate gamma-ray pulses from the neutron 
counter. The circuit was adjusted so that only neutrons over 800 kev were 
counted by the detector. This permitted the gamma rays from F'® to be 


*Obtained from Nuclear Enterprises, Winnipeg, Manitoba, Canada. 
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identified. For these measurements the neutron detector was placed at 0° to 
the incident Oj{ beam where, because of the large center of mass velocity, 
it could intercept a larger fraction of the reaction neutrons. 


EXPERIMENTAL RESULTS 


A portion of the gamma-ray pulse spectrum from the bombardment of H$ 
by a beam of 14-Mev O(§ is shown in Fig. 1. The high energy gamma rays 
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Fic. 1. Part of the gamma-ray pulse spectrum from the bombardment of H* by 14.0-Mev 
O}5 ions. 


from the H*(O}},ay)N® reaction are easily identified because this reaction 
has the highest Q-value (7.70 Mev) of the possible reactions emitting heavy 
particles. Lower energy gamma rays from the H*(Oj§,py)O'* (Q-value 3.73 
Mev) and H*(O/§,ny)F!8 (Q-value 1.28 Mev) are also observed. The gamma 
rays from F'’ were identified by measuring them in coincidence with neutrons. 
The gamma-ray pulse spectrum in coincidence with these neutrons is shown 
in Fig. 2. The well-known (Kuehner et al. 1958) 0.660, 0.94, 1.08-1.04 com- 
plex, 1.70, 2.10, and 2.53 Mev gamma rays can be observed in this spectrum. 





280 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


0-5 0 ek) 2:0 2:5 3-0 


DIRECT GAMMA-RAY PULSE SPECTRUM 


1000 


100 


COUNTS PER CHANNEL 


GAMMA-RAY PULSE SPECTRUM IN 
COINCIDENCE WITH NEUTRONS 





10 20. +30 40 50 6070 #80 90. 100 
CHANNEL NUMBER 


Fic. 2. A portion of the gamma-ray pulse spectrum from the bombardment of H® by 
16-Mev 03%, ions together with the spectrum in coincidence with fast neutrons from the reaction 
H*(O1$,,7)F'!8, The 1.29-Mev gamma ray is from the decay of A‘! produced by the Chalk River 
NRX reactor. This gamma ray is present during certain meteorological conditions. 


Evidence for higher energy gamma rays can also be seen. As a result of this 
measurement, which was carried out with 16-Mev Oj;{ ions, the incident 
oxygen ion energy was reduced to 14 Mev to remove the possibility of coinci- 
dences between gamma rays from the 3.725- and 3.790-Mev states in F'® 
interfering with the measurements on the gamma rays from O!’. In addition 
the Oj/§ beam was used for all subsequent measurements. 

To study the gamma rays from O!8 it was necessary to use the two large 
Nal crystals in coincidence. Gamma-ray coincidences from F!8 can be rejected 
because they lie in the lower channels of the 100-channel analyzer and in 
addition there are no known cascades in N” which could interfere with the 
measurements. Coincidence gamma-ray pulse spectra are shown in Fig. 3 
together with the relevant portion of the O!8 energy level diagram. In Fig. 3(a), 
the spectrum of gamma-ray pulses in coincidence with a gate set to include 
the total absorption peaks of the 1.98- and 1.94-Mev gamma rays is shown. 
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Fic. 3. Typical gamma-ray pulse spectra taken in coincidence. (a) The gamma-ray pulse 
spectrum in coincidence with a voltage gate set to include the 1.98- and 1.94-Mev gamma rays. 

(b) The gamma-ray pulse spectrum in coincidence with a gate set to include the 1.65-Mev 
gamma rays. The origin of these gamma rays is shown in the energy level diagram which 
includes the spin assignments made in this paper together with the parity assignments of 
Jaffe et al. 1960. 


Both crystals were at 90° to the beam and on opposite sides of the target for 
this measurement. 

In addition to the spectra shown in Fig. 3, spectra were taken with a voltage 
gate set above the 2-Mev peak. The resulting spectra showed that the 1.65- 
and 1.98-Mev gamma rays were not in coincidence with gamma rays greater 
than 2 Mev. Consequently the 1.65- and 1.98-Mev gamma rays cannot come 
from N, The pulses greater than 2 Mev in Fig. 3 require some comment as 
a calculation showed that they could not be random coincidences. It will be 
noticed from Fig. 2 that the voltage gate includes, in addition to the 1.98-Mev 
gamma ray, the 2.10-Mev gamma ray from F'8, As the 2.10-Mev gamma ray 
is in coincidence with neutrons there can be genuine neutron — gamma-ray 
coincidences between the two crystals. However, fast neutrons produce a 
continuum of pulses rather than a peak in the pulse spectrum and this is 
observed. The large velocity of the center of mass in the reaction causes the 
neutrons to be strongly peaked forward in the laboratory. Consequently when 
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the Nal(T1) crystal detecting the gamma rays at various angles to the incident 
beam is placed at 0° the pulses due to fast neutrons in coincidence with gamma 
rays, should increase in number. This effect was observed and it is shown 
for example in Figs. 5(a) and (0). 

To obtain the intensity of the gamma rays at each angle it was necessary 
to subtract the background due to fast neutrons. This was done by assuming 
that the spectrum of pulses induced by neutrons was the same above and 
below the gamma-ray peaks shown in Fig. 3. The subtraction of the tail of 
the 1.98-Mev gamma-ray spectrum line from the 1.65-Mev gamma ray was 
carried out using the known spectrum line shape of the 1.85-Mev gamma ray 
from Y*. 

Gamma rays of approximately 1.98 and 1.65 Mev appear in the coincidence 
spectrum. These gamma rays can be interpreted as resulting from cascades 
from the 3.92- and 3.63-Mev states in O'*. This interpretation is supported 
by the spectrum shown in Fig. 3(6). This shows the pulse spectrum from one 
crystal in coincidence with a gate set to include the total absorption peak 
of the 1.65-Mev gamma observed in the other crystal. The spectrum confirms 
the conclusions drawn from Fig. 3(a). As will be discussed below, there is a 
small contribution from the gamma rays from the 3.55-Mev state in O}8 
which results in coincidences between the 1.98-Mev gamma ray and 1.57-Mev 
gamma rays being observed at angles where the intensity of the 1.65-Mev 
gamma ray is low. 

It is important to consider the energy measurements made with the gamma- 
ray detectors because it is necessary to be sure that the 3.55-Mev state is not 
fed strongly at the bombarding energy used. The results of the measurements 
made by Silbert and Jarmie (1959) on the energy of the protons from the 
reaction O'*(H’,p)O'8 are given in Table I together with other measurements 
quoted by Jaffe et al. (1960). Jaffe et al. (1960) give a value of 82+7 kev for 
the energy separation of the 3.55- and 3.63-Mev states. Considerable difficulty 
was encountered in the measurement of the energies of the gamma rays in 
this experiment due to photomultiplier gain shifts which were a function of 
counting rate. Consequently, a counting rate meter was used to ensure com- 
parable counting rates when the reaction was being studied and when the 
calibration sources were used. Sources of Y** and Pr'** were used to provide 
1.850+0.008 and 2.185+0.015 Mev calibration gamma rays. 

These measurements can be analyzed in two ways: either by evaluating the 
energy separation of the peaks shown in Fig. 3(a), or by evaluating the actual 
energies. The first procedure is the more accurate and reliable since the energy 
separation is less sensitive to gain shifts resulting from changing counting rates. 
The measured energy separation of the two peaks was found to be 307+25 kev. 
The two calibration sources were used simultaneously to obtain the kev per 
channel. The energies shown in Table I can be used to calculate the energy 
separation of the cascade gamma rays from the 3.92-Mev state and those from 
the 3.55- and 3.63-Mev states. These calculated energy separations are 384+10 
kev and 302+10 kev respectively which imply that the 3.63-Mev state is 
excited more strongly than the 3.55-Mev state as the measured value was 
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307+25 kev. As will be discussed later there is evidence that the 3.55-Mev 
state is excited weakly in the reaction. 

An analysis of the width of the 1.98-Mev complex shown in Fig. 3(a) 
indicates that the energy separation of the two gamma rays is 20+13 kev. 
This latter energy separation is in approximate agreement with the value 
of 43+10 kev deduced from the energies shown in Table I. It should be pointed 
out that such an analysis is possible in this case because both gamma rays 
are equal in intensity. This is due to the inclusion of both gamma rays in the 
voltage window set on one of the counters at 90° to the beam and to the 
observation of the pulse spectrum shown in Fig. 3(a) at 90° to the beam in 
the other counter. 

The energies of the prominent peaks in the spectrum of Fig. 3(a) were also 
determined relative to the 1.850+0.008 Mev and 2.185+0.015 Mev calibra- 
tion gamma rays. These energies were found to be 1.677+0.025 Mev and 1.983 
+0.025 Mev. Because the upper peak is composed of two gamma rays which 
are close together in energy and which cascade from the 3.915+0.005 Mev 
state the expected peak position is 1.957+0.007 Mev which compares favorably 
with the measured value. The expected energy of the cascade gamma ray 
from the 3.634-Mev state to the 1.979-Mev state is 1.655+0.007 Mev which 
again compares favorably with the measured value. The errors on the measured 
energies do not include the possibility of an error due to the small difference 
in average counting rate between the reaction gamma rays and the calibration 
gamma rays. The fact that the measured energies are both about 25 kev 
higher than the expected values may be a result of this effect. 

The measurements discussed above were considered to be sufficient evidence 
that the 3.63- and 3.92-Mev states are preferentially excited at the incident 
energy and target conditions used. 

The crossover transitions from the 3.63- and 3.92-Mev states in O'’ were 
not observed in spectra such as that shown in Fig. 1. It could be estimated 
from such spectra that the crossover transitions were less than approximately 
15% of the cascade transitions. These crossover transitions could be searched 
for more easily by observing the gamma rays in coincidence with the appro- 
priate proton group. 

The angular correlations of the two complexes of gamma rays were measured 
in the geometrical arrangements shown in Figs. 4 and 6. The protons leading 
to the excited states of O'!8 were not observed. As was discussed above, the 
1.98-Mev complex consists predominantly of a 1.65-Mev gamma ray and a 
weak 1.57-Mev gamma ray. For the measurement of the angular correlations 
a voltage gate was set on the combined total absorption peaks of the 1.98- 
and 1.94-Mev gamma rays from one counter fixed at 90° and the coincident 
gamma-ray pulse spectrum observed at various angles (6,6) of the other 
counter. The correlations shown in Figs. 4 and 6 are corrected for the varying 
absorption of the target backing and target chamber walls which is a function 
of the angles (6,6). They were also corrected for a small effect due to the beam 
spot on the target not being at the center of rotation of the counter assembly. 
These corrections were estimated by observing the angular correlations of the 
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2.185-Mev gamma rays from a Pr! source placed in the target holder. Small 
corrections were also made for the aberration effect due to the motion of the 
source of gamma rays from the reactions. These corrections were in nearly 
all cases smaller than the statistical errors. The corrected angular correlations 
were fitted to a function of the form ao+d2P2(cos 6)+a4P.4(cos 6) or, in the 
case of Fig. 6(0), of the form 9+, cos 26+, cos 4¢. The coefficients ao, ae, 
a4, bo, bz, and 64 have in turn been fitted, as discussed in the next section, to 
the spin assignments 0 — 2 — 0 for the 1.65 Mev-1.98 Mev correlation and 
2— 2-0 for the 1.95 Mev-1.98 Mev correlation. The latter fit was made 
using a Datatron computer program written by one of the authors (A.J.F.). In 
measuring the correlations of Figs. 4 and 6, no attempt was made to establish 
the relative normalizations of the different distributions. The theoretical 
distributions, which are shown as solid lines, are normalized so as to be 
continuous at the junctions between them, and the same normalization factors 
have been applied to the experimental points. The results of the analysis are 
shown in Table II and Table III. 


TABLE II 


The experimental coefficients for the angular 

correlation of the 1.65- and 1.98-Mev gamma 

rays are compared with theory. For correlation 

(a) the coefficients refer to a correlation which 
is a function of an angle 90° — 4, 


Experiment Theory 


(a) a2/ao .28+ .13 .309 
a4/ao .59+.12 643 
(b) a2/ao 40+ .04 309 
a4/ao 61+ .05 643 
(c) a2/ao .214+.06 0 
a4/do .03+ .07 0 


TABLE III 


The experimental coefficients for the angular correlation of 

the 1.98- and 1.94-Mev gamma rays are compared with 

theory for the quadrupole—dipole amplitude mixture ratio 

x = 0.1. The experimental correlations are shown in 
Fig. 6 together with the theoretical curve 





Experiment Theory 

(a) d2/ao —0.17+.03 —0.158 
a4/ao 0.10+ .03 0.105 

(b) be/bo 0.40+ .03 .424 
b4/bo 0.06+ .03 .002 

(c) a2/ao 0.58+ .06 .529 
a4/do 0.12+ .08 .002 


ANALYSIS OF THE EXPERIMENTAL DATA 


(a) The Correlations of the Gamma Rays from the 3.63-Mev State 
The strong angular correlations shown in Figs. 4(a@) and (6) immediately 
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suggest the assignment of spin 0 to the 3.63-Mev state in O'8% The solid 
curves shown in Fig. 4 are theoretical curves, for the spin sequence, 0 — 2 > 0, 
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Fic. 4, Gamma-ray angular correlations taken with three geometrical arrangements of the 
counters shown in the figure. The energies of the gamma rays observed in each counter are also 
shown. The solid lines on the angular correlations (a) and (c) are for the spin sequence 0 — 2 


=> @, 


attenuated for the finite sizes of counters employed. Since the initial state of 
the cascade is assumed to have spin zero, the gamma-ray correlation is that 
of a decaying radioactive nucleus, and is given by the equation 
(1) W = 1+(Q2)? 0.357 P2(cos $1) + (Qs)? 1.143. P4(cos ¢1). 


The axis of quantization is defined here by the fixed counter and the target, 
and is different from that assumed below in the analysis of the correlations 
from the 3.92-Mev state. ¢; of equation (1) is defined as shown in Fig. 4(d). 
The same equation describes the correlation of Fig. 4(a@) with ¢; replaced by 
90°—6,. The attenuation coefficients Q. and Q4 were taken from the tables 
of Gove and Rutledge (1958) and Rutledge (1959). The values used in Figs. 4(a) 
and (c) were Q2 = 0.93 and Q, = 0.75. The agreement is very good in Figs. 
4(a) and 4(b). However, there is disagreement in Fig. 4(c) because there the 
predicted correlation is isotropic. The reason for this discrepancy was apparent 
from the gamma-ray spectra taken for the correlation shown in Fig. 4(c). Two 
gamma-ray pulse spectra for this correlation are shown in Figs. 5(a) and 5(0). 
At 6, = 0° in Fig. 5(6) the 1.65-Mev gamma-ray total absorption peak is 
broader than that recorded at 0; = 75° which is shown in Fig. 5(a). The 
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Fic. 5. Typical pulse spectra taken during the angular correlation measurements. The 
voltage gate was set to include the 1.94- and 1.98-Mev gamma rays. For the spectra one counter 
was situated at @ = 90°, ¢ = 180°, and the other was moved in the ¢ = 90° plane. This angular 
correlation is shown in Fig. 6. These spectra show the broadening of the 1.98-Mev complex 
due to the 1.98- and 1.94-Mev gamma rays from the 3.92-Mev state. In spectrum (6) the 
1.65-Mev complex is perceptably broadened over that in (a) indicating the presence of a weak 
1.57-Mev gamma ray from the 3.55-Mev state in O#, 


broadening from 6.5 channels in 5(@) to 8.5 channels in 5(6) can be explained 
by the appearance of the 1.57-Mev gamma ray from the de-excitation of 
the 3.55-Mev state in O'8 with an angular correlation such that it is much 
stronger at 0° than at 75°. A rough estimate of 25% for the intensity of the 
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1.57-Mev gamma ray relative to the 1.65-Mev gamma ray at 6; = 0° was 
obtained from measured spectra. Admixtures of this magnitude of a 4 - 2 > 0 
correlation having the lower magnetic substates of the 4+ state the most 
strongly populated are compatible with the anisotropy appearing in Fig. 4(c). 
The effect of such admixtures on the correlations of Figs. 4(a) and 4(d) will 
not be great enough to change the conclusions. The good agreement between 
theory and experiment shown in Fig. 4(@) and Fig. 4(b) together with agree- 
ment shown in Table II can be taken as strong evidence for the spin 0 assign- 
ment for the 3.63-Mev state in O'8. In Table II the finite value of the corre- 
lation coefficient listed in Section (c) is attributed to the presence of the weak 
1.57-Mev gamma ray as discussed above. 

A search for fits to the data assuming the spin sequences 1 > 2-0, 
2-2-0, 3-2-0, and 4 > 2 — 0* was made with the aid of the com- 
puter program described briefly below. The theoretical Legendre polynomial 
coefficients given by the best of these were grossly discordant with the 
measured values and it is concluded that the measured correlations exclude 
these assignments. 


(b) The Angular Correlations of the Gamma Rays from the 3.92-Mev State 

The angular correlations of the 1.98+1.94 Mev gamma rays are shown in 
Figs. 6(a), (6), and (c). The solid curve shown is the theory for the spin 
sequence 2— 2-0 with the first radiation pure dipole. The theoretical 
angular correlations were calculated in this case using the procedure outlined 
by Batchelor et a/. (1960) and discussed in some detail by Litherland and 
Ferguson (1960). As shown by Litherland and Ferguson (1960) the angular 
correlation of the gamma rays from an axially symmetrical state J;, is 





(2) W wv (01,02,0) = ood (—1)***70 OT (8) x? Dy ay (tt) Xx a4 (1,02, 6). 





In equation (2) the state of spin J; decays to a state of spin J, by emitting 
2¥1-pole radiation. The state of spin Jz decays to a state of spin J by emitting 
2¥2.pole radiation. The first gamma ray is observed at a polar angle 6; and 
the second at 6. The azimuthal angle between the two gamma rays is ¢. 
The polar axis, i.e. the axis of quantization, lies in the beam direction. The 
coefficients Dx(tt’) have been tabulated by Ferguson and Rutledge (1957). 
The functions Xx%17(61,02,6) are defined by Ferguson and Rutledge (1957). The 
attenuation coefficients Qx, Qy have been tabulated by Gove and Rutledge 
(1958) and by Rutledge (1959). The quantities 7(s) are used to specify the 
alignment of the state J; and are related to the populations of the magnetic 
substates of the state J; by the equation 
*It can be shown that the three measured ‘ ‘geometries’ ’ provide only six independent constants. 
Fits made assuming the spin sequences 3 — 2 — 0 and 4 + 2 — 0 involve seven parameters. 
For the 3-+2-0 case these consist of one E2/M1 ratio, three population parameters, 
P(m) (see below), and three normalization factors. For the 4 2—0 0 case no E2/M1 ratic 
is present, but there is an additional population parameter. Thus unique fits cannot be expected 
for these cases. However, two normalizations may be determined at least approximately, for 


example, by joining the curves at the common points. To attempt the3 — 2 — Oand4 — 2 — 0 
fits mentioned, the normalizations of Figs. 4 and 6 were used. 
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Fic. 6. Gamma-ray angular correlations taken with the three geometrical arrangements of 
counters are shown in the figure. The two gamma rays of the 1.98-Mev complex could not be 
resolved so that the correlations of their sum was measured. The solid lines on the angular 
correlations are for the given sequence 2 — 2 > 0. 


(3) P(m) = ie T(s)(sJim—ml|l0)’. 


x is the E2/M1 amplitude ratio of the first gamma ray and p is an exponent 
having the value 0, 1, and 2 for the dipole-dipole, dipole—quadrupole, and 
quadrupole-quadrupole terms respectively.* In equation (3) the quantum 
numbers satisfy the relation s+1 = J, where the angular momentum quantum 
number / is chosen so that the Dxs,(tt’) coefficients can be obtained from the 
tables of Ferguson and Rutledge (1957). For the calculations discussed in 
this paper it was found convenient to choose / = J,. 

In the correlations shown in Fig. 6 one of the counters was fixed at 90° to 
the beam direction. For the arrangement of counters shown in Figs. 6(a) and 
(c) equation (2) simplifies to 


(4) Wir (0) = LL T(s)x?A(LL)P, (cos 8) 


where 


AAEL) = 2 ore Ox QuDrar (tt’). 


*Since the second gamma-ray transition involves the spin sequence 2 — 0, this transition 
is a pure E2 and no mixing parameter is required for it. 
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The coefficients atxy have also been tabulated by Ferguson and Rutledge 
(1957). For the case shown in Fig. 6(6) equation (2) simplifies to 


(5) Wi (o) = 2) Dd T(s)x’By(LL’) cos No 


where 


By(LL’) = 2, ox Ox OuDieae(t'). 
The coefficients 4,(LL’) and By(LL’) are tabulated in Table IV of this paper 


TABLE IV 


The coefficients A,(11) and B,(11) of the angular correla- 

tions are tabulated for the sequence of spins 2 + 2 — 0. 

The coefficients contain the attenuation factors Qxr 
discussed in the text 


s 0 1 2 


(a) Ao(11) 2.5325 1.8161 2.3825 
A(11) —0.4135 0.2905 —0.7319 
A,(11) 1.2594 —0.6407 0.2746 


(b) Bo(11) 0.8912 0.6740 1.2157 
B11) 0.7146 0.0413 0.2101 


(c) Ao(11) 1.3231 1.3578 1.9715 
A,(11) 1.9893 0.1521 —0.0607 
A,(11) 0.0661 —0.0440 0.0189 


for the case L = L’ = 1 and for J; = 2, Je m2, I = 0, and } = 2. The 
coefficients apply only to the case where the two gamma rays L, and Lz are 
not distinguished by the gamma-ray detectors because of their being so close 
in energy. 

Preliminary calculations revealed that the measured Legendre polynomial 
coefficients of Table III are compatible with the spin sequence 2-2-0 
and with a small admixture of quadrupole radiation in the 1.94-Mev tran- 
sition. This work has been refined by a least squares fit carried out with the 
Chalk River Datatron, the resulting coefficients being listed in the ‘‘theoretical’”’ 
column of Table III.* The adjustable parameters of this fit were the ratios 
T(1)/T(O), T(2)/7(0), the quadrupole-dipole amplitude ratio, and, since the 
relative normalizations of the three different cases were not measured, a 
normalization constant for each of them. These constitute six parameters to 
be determined by the nine measured Legendre polynomial coefficients: In 
Fig. 6 the three theoretical correlations have been renormalized so as to be 
continuous at the junctions between them and the same normalization con- 
stants have been applied to the experimental points. 

The results of the fit are as follows: the E2/M1 amplitude ratio of the 

*It will be noted that two least squares fittings have been made in succession here, the first 
to obtain the ‘‘measured” Legendre polynomial coefficients, the second to obtain the “‘theo- 
retical’’ ones. While there are some minor objections to this double process, it is more con- 


venient to program on the computer and more economical of computer time than a program 
which would fit directly the measured data in terms of theoretical parameters. 
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1.94-Mev transition is 0.1.1, the relative populations of the magnetic 
substate of the initial 2+ state, P(Q): P(1): P(2) = 1:0.67+0.17:0.08+0.07. 
The presence of alignment here implies that the orbital momentum of the 
incoming particles is greater than zero, and that the compound state has 
a spin greater than 1/2. The best tits available assuming the spin sequences 
1-2-0, 3-2-0, and 4—2->0, were also obtained. For these the 
theoretical Legendre polynomial coefficients differed from the measured ones 
by considerably more than the experimental error, so that these possible 
assignments are eliminated. The sequence 0 > 2-0 is eliminated because 
the correlation of Fig. 6(c) is not isotropic. 


CONCLUSIONS 


The experiment described above has obtained further information on the 
properties of the excited states of the O'§ nucleus at 3.63 Mev and 3.92 Mev. 
In addition the measurements reported by Gove and Litherland (1959) on 
the properties of the 3.55-Mev state are unaffected by the presence of the 
then unknown 3.63-Mev state. This point is discussed in detail in the Appendix I 
to this paper. 

Measurements of the angular correlations of the cascade gamma rays from 
the 3.63-Mev and 3.92-Mev states in O'8 have permitted the assignment of 
spins 0 and 2 respectively to those states. Also the quadrupole amplitude 
divided by the dipole amplitude, x, of the 1.94-Mev gamma ray de-exciting 
the 3.92-Mev state was estimated to be 0.140.1. The information obtained 
on O'8 as a result of this experiment together with information previously 
obtained is summarized in Fig. 7. The positive parities of the 3.63- and 3.92- 
Mev states have been inserted because of the recent evidence from O'*(t,p)O18 
double stripping distributions obtained by Jaffe et al. (1960). 

The ratios of the energies, and the values of the spins and parities of the 
first four states in O'8 are remarkably similar to those found for the vibrational- 
like states found in medium weight nuclei such as Cd!“ (Motz 1957). However, 
the well-studied vibrational states in medium weight nuclei exhibit enhanced 
electric quadrupole transitions which have not yet been found in O'8, with 
the possible exception of the 4+ to 2+ pure £2 transition between the 7.13- 
Mev and the 1.98-Mev levels which shows some enhancement (Gove and 
Litherland 1959). The Weisskopf estimate (Wilkinson 1956) of the ratio of 
the electric quadrupole width to the magnetic dipole width is ~10- for the 
1.94-Mev gamma ray de-exciting the 3.92-Mev state in O18. However, if O'8 
is considered to be two neutrons attached to the closed shell nucleus O", 
electric quadrupole radiation can arise only from the recoil of O'* during the 
gamma-ray transition. The Weisskopf estimate is then multiplied by Z?/A‘ 
(Wilkinson 1956) which lowers the ratio to ~6X10~7. The experimental 
amplitude ratio of +0.1+0.1 from the angular correlations shown in Fig. 6 
corresponds to an intensity ratio of <.04 which implies an enhancement 
of <10°. It is interesting to recall that the 872-kev E2 transition from the 
first excited state of O" is strongly enhanced over the Weisskopf estimate 
multiplied by Z?/A‘*. 
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Fic. 7. The known energy levels of O'§ with the assignments and gamma rays studied in 
this experiment included. 











In conclusion a more accurate experimental determination of the dipole- 
quadrupole mixing ratio of the 1.94-Mev gamma ray would be of considerable 
interest together with a determination of the lifetimes of the 1.98-Mev, 
3.63-Mev, and 3.92-Mev states in O'’. A careful search for the de-excitation 
of the 3.92-Mev states by means of a ground state transition would also be 







of value. 


* 
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APPENDIX I 
SPIN ASSIGNMENTS FROM C!4(qa,y)O!8 


In the C'(a,7)O!® reaction at the 1.14-Mev resonance it was assumed by 
Gove and Litherland (1959) that there were only two primary transitions, 
that to the level at 1.98 Mev and that to the 3.55-Mev level. With the present 
information that there is a spin-0 level at 3.63 Mev the possibility exists 
that both it and the 3.55-Mev level are fed by a primary transition from the 
resonance in question and that these two primaries of energy 3.48 and 3.58 
Mev respectively would be unresolved in the direct spectrum. Coincidence 
measurements at this 1.14-Mev resonance proved that the level or levels 
near 3.6 Mev fed by primaries do not de-excite directly to ground thus avoiding 
the possibility that more than two gamma rays are unresolved at energies 
near 3.6 Mev in the direct spectrum. 

It was shown by Gove and Litherland (1959) in an appendix that the 
angular distribution of the primary transition from the 1.14-Mev resonance 
to the 1.98-Mev first excited state requires that the resonance has J = 2+ 
or 4+. If J = 4+ there is no problem since there will be no primary tran- 
sition to the spin-0 state at 3.63 Mev. However, if the resonance has J = 2+ 
the situation is more complicated. It is easy to show that the resonance 
cannot be 2+ by considering what this would imply about the direct angular 
correlation of the 1.98-Mev gamma ray. As before (Gove and Litherland 
1959, Appendix) the correlation of this gamma ray can be corrected for the 
part due to direct feeding of the 1.98-Mev state by a primary transition from 
the 1.14-Mev resonance. When the latter has J = 2+ the corrected corre- 
lation is W(@) = 1+(0.90+0.15) P2—(0.50+0.21) Py. One can assume that 
the 3.55-Mev state has any spin value between 0 and 4. Higher spin for the 
state would result in a vanishingly low probability for a transition from a 
J = 2+ resonance. Spin 0 can be eliminated since a correlation of the type 
2+ (£2) 0+ (£2) 2+ (£2) 0+ with the first two gamma rays unobserved is 
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spherically symmetric. For spin 1 there would be no P, term in the correlation 
whereas a large one is required experimentally. For spin 2, one can show that 
the maximum value of the coefficient of P2, is +0.178 for any combination of 
dipole, quadrupole mixtures in the first two radiations. This is already too 
small and any feeding of the spin-0 state at 3.63 Mev adds in a spherically 
symmetric distribution which makes it even smaller. Similarly for spin 3 the 
maximum P? coefficient is +0.49 which is again too small. There remains only 
the possibility that the 3.55-Mev level has spin 4 with spin 2 for the 1.14- 
Mev resonance. In this case the theoretical correlation 2+ (£2) 4+ (£2) 
2+ (£2) 0+ with the first two gamma rays unobserved is W(6) = 1+0.40P, 
—0.14P,. If the 3.63-Mev state is fed as well one must add a spherically 
symmetric distribution to this which reduces both coefficients and puts them 
even more outside the required values. If the 3.63-Mev state is not fed, the 
arguments employed by Gove and Litherland (1959) do not permit the com- 
bination spin 2 for the 1.14-Mev resonance and spin 4 for the 3.55-Mev level. 
Hence the spin assignments made by Gove and Litherland (1959) are unchanged 
by the discovery of a spin-0 level at 3.63 Mev. 


APPENDIX II 


Although the triton energy in the center of mass system is only about 2.2 
Mev, when 14-Mev O/$ ions bombard H*, the center of mass has a velocity 
of about 1.1 10° cm per second in the laboratory. This implies large (~3%) 
doppler shifts of the gamma rays from the residual O'8 nuclei because they 


will be travelling at approximately the velocity of the center of mass. The 
large velocity of the recoiling nuclei from reactions with light elements induced 
by O' ions makes it possible in principle to use a technique described by 
Devons, Manning, and Bunbury (1955) to measure the lifetimes of the resulting 
excited nuclear states. They showed that it is possible to estimate the lifetime 
of an excited state produced in a nuclear reaction by observing the attenuation 
of the doppler shift when the recoiling nucleus is slowed down in a solid. The 
method is most successful when the slowing down time in the solid is of the 
order of the lifetime to be measured. For example, using the data given by 
Manning (1955) an O" ion travelling initially at 10° cm/sec in copper comes 
to rest in approximately 10~" sec. In the experiment described in this 
paper this method was not exploited because suitable targets were unavailable. 
It was possible, however, to obtain a lower limit for the lifetime of the 
3.63-Mev state in O'8 from the measured spectra. 

The gamma-ray spectra which were used to determine the angular distri- 
butions shown in 4(6) and 4(c) were scrutinized for energy shifts. No significant, 
<1%, shift for either gamma ray was observed for the sequence of spectra 
taken to obtain Fig. 4(b). This is to be expected as both gamma-ray counters 
are at 90° to the incident oxygen beam. The spectra upon which Fig. 4(c) is 
based were, however, taken with one of the gamma-ray counters at a variable 
angle 6 to the beam. Again no significant, <1%, shift between 6 = 0° and 
6 = 90° was observed for the 1.65-Mev gamma ray but a shift of ~1.8+1.0% 
was observed for the 1.98-Mev complex of gamma rays. 
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As the target was a thick one it is not possible to estimate the average 
doppler shift before the recoiling O'* ions start to slow down in the backing. 
However, because the Coulomb barrier to the H*+O!* reaction causes th. 
yield of the reaction to fall rapidly as the energy of the oxygen ions is decreased 
the average doppler shift is expected to be near the maximum value of about 
3%. Consequently the <1% doppler shift observed for the 1.67-Mev gamma 
ray implies that most of the recoiling nuclei have slowed down before the 
gamma rays from the 3.63-Mev state in O'8 are emitted. The slowing down 
time for an O} ion of initial velocity 10° cm/sec in copper is approximately 
10-2 sec which is therefore an approximate lower limit to the lifetime of the 
3.63-Mev state in O'%, 

The ~1.8+1.0% shift of the 1.98-Mev complex of gamma rays is con- 
sistent with the 1.94-Mev M1 gamma ray showing the full doppler shift of 
~3% and the 1.98-Mev £2 gamma ray showing no doppler shift. 

It should be emphasized that these estimates are included merely as a by- 
product of the data taken for the angular correlations. Much better estimates 
of the lifetimes could be obtained if it was possible to obtain thin tritium 
targets (a) in zirconium thin to the incident oxygen ions, (6) in thin zirconium 
backed with calcium, and (c) in thin zirconium backed with copper. 





AN ITERATIVE METHOD FOR NEUTRON TRANSPORT 
PROBLEMS WITH SPHERICAL SYMMETRY! 


W. R. ConkKIE? 


ABSTRACT 


An iterative method developed previously for transport problems with plane 
symmetry has been extended to problems with spherical symmetry. Particular 
application has been made to the problem of the neutron distribution for a black 
sphere embedded in a purely scattering medium with sources at infinity. The 
results compare favorably with those of other workers for this problem. 


1. INTRODUCTION 


For many of the problems of interest to reactor physicists the stationary 
neutron distribution in a reactor can be mathematically described as being a 
solution of the transport equation (Davison 1957) 

(1.1) 7. WN (roa) + EO) = Jw W(ro'a!ydv'aa’ + S(r,vQ). 
tot 

In this equation V(r,vQ) is the probable number of neutrons in the volume 
element dV around the point r, with velocity V = vQ, where v is the magnitude 
of the velocity. /,.,(v) is the total mean free path, and K(v,v’) is the differential 
scattering probability for scattering from v’ to v. S represents any sources 
which are present and we have assumed cross sections independent of r. 

It is usually very difficult to obtain accurate solutions to this equation and 
it is necessary to turn to approximate methods of solution which in many cases 
have only a limited range of application. For example, for slowing down 
problems the Greuling—Goertzel method (Greuling e¢ a/. 1953) has been 
useful, and for thermalization problems a method combining the spherical 
harmonics method with polynomial approximations in the velocity variable 
has been applied to a simple problem with plane symmetry (Conkie 1960a). 
However, these methods and others suffer from their restriction so far to a 
P; approximation (Davison 1957, Chapter X ef seg.). That is, the neutron 
distribution consists only of a spherically symmetric term plus a term linear 
in the directional cosines of the angular variable 2. Improvements in this 
approximation can in principle be made. For example, for the thermalization 
problem mentioned above (Conkie 1960a) the extension to higher Py approxi- 
mations is straightforward but will involve a large amount of computational 
labor. Also, for many problems the usual form of the spherical harmonics 
method is slowly convergent and this may not be the best way to improve the 
accuracy of these various approximate methods. 

Since it is relatively easy to find low order spherical harmonics approxima- 
tions for various problems we might consider iterative methods of improving 
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these solutions, since the operators forming the transport equation (1.1) are 
easy to apply to the solutions usually obtained. 

To this end of developing more accurate solutions of the transport equation 
from low order spherical harmonics approximations, the author has been 
developing an iterative method for solving the transport equation, and has 
already applied this method to the Milne problem of one velocity-group trans- 
port theory or radiative transfer (Conkie 1959). 

The present work still remains within the limits of one velocity-group theory, 
but extends the iterative method to problems with spherical symmetry. We 
will consider in particular the problem of the neutron distribution around a 
black sphere embedded in a purely scattering medium with neutron sources 
far from the sphere. This problem has been treated before by various methods 
(Davison and Kushneriuk 1946; Davison 1951) including the spherical 
harmonics method. This method does not give accurate results for a sphere of 
small radius even in high order approximations. We will focus attention on the 
calculation of the linear extrapolation length, A, which is defined as 


Pas (r) 

cee » = W/dr)oma 7) | rae 

where p,s(7) is the asymptotic neutron flux and a is the radius of the sphere. 
The parameter d is useful in obtaining the correct ratio of flux to its derivative 
in the diffusion theory approximation. It should be noted that although we 
will calculate only in detail, the complete neutron distribution should be 
obtained from our results with an accuracy equal to that with which we obtain 
the values of the extrapolation length. 

In Section 2 we recapitulate the main features of the iterative method applied 
to problems with plane symmetry. This will serve to illustrate the concepts 
necessary for carrying out the procedure, and the approach will be somewhat 
different from the previous exposition (Conkie 1959). 

In Section 3 we extend the method to a problem with spherical symmetry, 
the black sphere problem, with emphasis on the linear extrapolation length. 

Section 4 contains the comparison of the results of Section 3 with those of 
other workers and the conclusions derived from these comparisons. 

In Section 5 we discuss the extension of the method to general multilayer 
problems with spherical symmetry. 

Section 6 contains conclusions as to the general applicability of the method. 


2. PROBLEMS WITH PLANE SYMMETRY 


For problems with plane symmetry the one group source-free transport 
equation with isotropic scattering is (Davison 1957) 


(2.1) pte) yen) =f veu' du. 


¥(x,u) is the angle-dependent flux, usually called the angular distribution, 
u is the cosine of the angle between the x direction and the direction of velocity, 
c is the mean number of secondaries per collision, and distances are measured 
in units of the mean free path. 
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The flux, p(x), is given by 


p(x) = ¥ ¥(x,u)du. 


It is possible to solve equation (2.1) by the spherical harmonics method, 
which for this problem consists of expanding the angular distribution in the 
form 


vem) = De vale) EYP) 


where the P, (x) are Legendre polynomials and the expansion is truncated 
at the Vth polyno 

By this means the integrodifferential equation (2.1) is replaced by a system 
of (V+1) coupled differential equations, which are readily solved for the 
moments (x) in the form 


Ya(x) = a &n(v;)e"™* 


where the g,(v,) and the exponents v,; depend on the order of approximation 
attempted (Davison 1957, Chapter X). 

It should be noted that there are two main sources of error in applying the 
usual form of the spherical harmonics method. First, there is the error of 
truncating the infinite system of differential equations to (V+1) equations. 
Secondly, instead of satisfying the exact boundary conditions of the problem 
being considered, with the complete neutron distribution being continuous 
across an interface, we use the approximation that only certain angular 
moments of the distribution be continuous. 

We can expect, in deriving solutions in this approximation, that the flux is 
determined more accurately than the higher angular moments, since it is an 
integral over the approximate angular distribution and presumably the error 
in the angular distribution is oscillatory because of our orthogonal functions 
expansion. 

This suggests putting the Py approximation to the flux, py(x), in place of 
the integral in equation (2.1) and simply solving the resulting differential 
equation. 


(2.2) pone 


t+ (xu) = 5 Px (2). 


This will give us an ue result for ¥(x,u), which can be made to satisfy 
the exact boundary conditions of the problem. This improvement of the 
distribution at the boundary is illustrated for the particular case of the Milne 
problem by Kourganoff (1952, p. 99). This type of iteration has been considered 
by various other authors previously and it was found in most cases that 
although ¥(x,) was improved considerably, the flux, p(x), was not improved 
very much in one iteration. That is, it was mostly the higher moments of the 
distribution which were improved. The main effect of this method seems to be 
that of improving the angular distribution, especially at the surface of the 













298 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 





region being investigated. In other words this scheme will largely compensate 
the second source of error mentioned above, but the errors due to truncation 
of the infinite system of equations still remains. This last statement is con- 
firmed by the fact that the exponents v,, which depend strongly on the order 
of approximation, remain unchanged by this iteration. 

Let us now consider a method of solving equation (2.1), which involves no 
truncation. To this end we will consider the Fourier transformed form of 


equation (2.1) 










(2.3) (1—ink) (em) =$ fo on! du+ Su) 










where 


ll 


olen) = [ev eu)dx 


S(u) = »¥(0,u)—vv(a,n), 


where for definiteness we have considered the neutron distribution in an 
infinite slab, bounded by the planes x = 0 and x =a. We have defined 
¥(x,u) to vanish outside the slab, and consequently the surface terms which 
make up S(u) appear as source terms in equation (2.3). This technique for 
replacing a finite medium problem by an infinite medium problem is discussed 
more fully by Case, deHoffman, and Placzek (1953, Section 17).* 

Equation (2.3) can now be solved for ¢o(&), the Fourier transform of the 
flux p(x), by dividing through by (1—iuk) and integrating over uw. We obtain 
then 














a Pe a 
(2.4) do(k) = yi k tan ky [. io du 







p(x = 5 J do(k)e” "dk. 
Tv “ca 





If S(u) is known, the integration in equation (2.4) can be done and the 
inverse Fourier transform of ¢o(&) taken using the methods of contour integra- 
tion. The type of contour used is shown in Fig. 1 (CHP, Appendix C). If the 
contour is taken in the upper half plane, and c < 1, [l1—(tan~'k)/k]-' has a 
pole at k = ixo and a branch point at k = 7. Thus we get a contribution 
from the pole at k = ixo and from integrating along both sides of the branch 
cut extending from k = 1 to k = i1~. The result for p(x) is obtained in the form 


p(é) a Pas (€) + Pnon-as (é), 


where p,s(&), the so-called asymptotic part of the solution is given by 












*Hereafter this reference will be called CHP. 
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i 
0 
Fic. 1. The form of contour required for evaluating p(x) from go(k) given by equation (2.4). 
OKo i S(u)du Kot 
2.5 = -—¢— a GF HOF 
(2.5) Pas (€) ac J_1 (=F xo) 


Pnon-as(£) is the non-asymptotic solution, sometimes known as the transient 
part of the solution, since it represents terms which die off rapidly as we move 
away from a source or boundary, and is given by 


(2.6) Pnon-as (£) = *ave(o)e ty SH) ; SW')— St) an’ 
0 M —y 6 eomp) 


In these equations the upper and lower signs hold for — > 0 and — < 0 respec- 
tively and = x or (x—a) depending on the term of our particular S(u) which 
is being used in the above integrals. g(c,u) is defined in CPH, Section 14, 
and is also given in Appendix ITI. 

We see then that as long as S(u) is known we can find our solution. Un- 
fortunately S(u) is not known, but depends on the solution of equation (2.1), 
which we are trying to find. However, since we know that solving equation 
(2.2) with accurate boundary conditions leads to much improved surface 
distributions, we can try using the improved S(u) calculated from equation 
(2.2), in equations (2.5) and (2.6). 

The iterative method which we propose consists then of three steps. 

(1) We determine a low order spherical harmonics approximation to the 
neutron flux. 

(2) We use this flux as a source term in equation (2.2), whence we improve 
the angular distribution, satisfying the boundary conditions of our problem 
accurately at the various surfaces. 

(3) We use these improved surface distributions as source terms in the 
Fourier transformed transport equation from which we can obtain an improved 
neutron flux. In practice this amounts to using the isotropic and anisotropic 
plane and point source distributions of CHP, together with the theorem given 
in Section 17 of that work. 

With respect to the two sources of error in the usual spherical harmonics 
method discussed previously, we see that step 2 attempts to improve the 
accuracy with which the boundary conditions are satisfied and thus to give 
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better surface terms. In step 3 the improved surface terms are then used in a 
complete solution of the transport equation which involves no truncation 
errors. 

It should be noted that by introducing step 3 we get the correct form of 
solution as given by equations (2.5) and (2.6) with the exact exponent, or 
reciprocal diffusion length, in the asymptotic solution and the continuous 
spectrum of exponents, 1/u, in the non-asymptotic solution. 

In principle steps 2 and 3 could be repeated to carry out: more iterations, 
but for most problems this would become laborious and we shall see that good 
accuracy can usually be obtained from just one stage of this process. 

As an example of this iterative process we will consider the calculation of 
the neutron distribution in the Milne problem. The results for the quantity 
called the extrapolated end point, 2, are easily calculated in various Py 
approximations of the standard spherical harmonics method. A description 
of the Milne problem and definition of the extrapolated end point together 
with various Py results for this quantity are given by Weinberg and Wigner 
(1958). The results for 29 for several Py approximations for a pure scatterer 
(c = 1) are given in Table I. 

























TABLE I 
Values of zo for various Py approximations. c = 1.0 













N 1 3 5 7 15 
Zo .5774 -6940 - 7039 . 7069 - 7096 










The exact value to four figures is 2) = .7104. 

The details of the iterative method applied to the P; solution have been 
given elsewhere (Conkie 1959) and the result is 29 = .7113. It is also found 
that other parameters of the distribution obtained from the iterative method 
are given with equal accuracy, which should be sufficient for most practical 
purposes. 

We see then, that by iterating the P; solution, which is easily obtained, we 
can obtain results which are as accurate as the P;s5 results for this problem. 
It should be appreciated that there is a great amount of computing involved 
in calculating the P;5 solution for this relatively simple energy-independent 


problem. 


















3. CALCULATION OF THE FLUX FOR THE BLACK SPHERE PROBLEM* 


Instead of developing the iterative method for general problems with spheri- 
cal symmetry, we will consider the specific problem of a black sphere embedded 
in a purely scattering medium with neutron sources far removed from the 
sphere. We will concentrate on the calculation of the linear extrapolation 
length, A, defined by equation (1.2), although, as will be seen, the complete 
neutron distribution can be obtained from our results. 

The P, result for the flux for this problem is readily calculated (Mark 1957), 










*A preliminary account of this work was presented at the June 1960 meeting of the American 
Nuclear Society (Conkie 19605). 
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and using Mark’s boundary conditions that the angular distribution at the 
surface of the sphere should vanish at the positive root of P2(u) = 0, we get 






(3.1) nie) = (442) -2, N=1 


a r 





where a is the radius of the sphere expressed in units of the mean free path of 
the scattering medium and gp = 1/+/3 = .5774. Then 


= [go/a*]/[1/a?] = go = .5774, for all a. 










Carrying on to step 2 of our iterative procedure, we will use py(r) given by 
equation (3.1) as a source term in 





(3.2) Q.Vyv(r,2Q)+yY(r, 2) = = five, Q’)dQ’, 





which leads to the equation 













(3.3) 2. V(r.) +¥(1,Q) = Fpw(r), 









for our problem with spherical symmetry. 
The solution of this equation, which we will call the improved angular 
distribution is given by 














V imp (F, Q) = ¥ imp (7,4) a et £ J, ton(e-Ra)ar, 





(3.4) 





where p =Q.1/r. 
In our case, with c = 1 





oie 
‘ce dR \ 
(3.5) V imp (7,4) = Aft +o . V(R?+r'—2rRu) 





with the surface condition 


(3.6) Vimp (ttt) [rae = 0, o> 









Since we will need only the results from equation (3.5) at the surface of the 
sphere, where condition (3.6) holds, an expansion of ¥(a,u) in shifted Legendre 
polynomials, adjusted to the range —1 < » < 0 (Lanczos 1956) is suitable for 
further calculations. 

In practice only the third term on the right-hand side of equation (3.5) 


needs to be expanded to give 
si e "dR 
0 V(R’+a°—2aRz) 


where the coefficients C,(a) are obtained from 









(3.7) wa > (2n+1)C,(a)P*(u), whips 









eR 
(3.8) Cala) = i Paluldu J, Ta? 2aRu) 
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These integrals were evaluated using high order Gauss—Legendre quadrature 
for the w integration and high order Gauss—Laguerre quadrature for the R 
integration. The computations were done on the Chalk River Datatron and 
results were obtained for various values of a. The rapid convergence of the 
expansion (3.7) which occurred for all cases considered is illustrated by the 
results shown in Table II for a = 1.0. 


TABLE II 
Coefficients C,(a) for a = 1.0 





‘ 


n 0 1 
C,(1) 6.64432 X 107 2.58298 X 107? 2.16066 X 10-* 


n 3 4 5 
Cr(1) 2.31098 X 10-4 2.76412 x 10-5 3.55770 X 10-6 





Let us now consider the Fourier transformed form of equation (3.3), with 
(3.9) ¢(k,Q) = f e**y(r. Q)dV, 


where the integral is over r space. By multiplying equation (3.3) by e** and 
integrating over 7, we obtain 


(3.10) (1-i@. k)$(k,Q) = ~ | $(k, 2’)do" 
dr 


_ fa .dSe** y(r.Q) 
where the last term, obtained by Gauss’ theorem, is a surface integral over 
the sphere surface. Q.dS = a?(Q.r/r)dQ,. 

Since ¥(r,Q) is a function only of 7 and (Q.r/r), ¢(K,Q) is a function only 
of k and (Q@.k/k), where & is the magnitude of the vector k. 

Calling the last integral on the right-hand side of equation (3.10), X (k,Q) 
= X(k,uo), where uo = Q .K/k, we obtain, by dividing through by (1—7Q. k) 
in equation (3.10) and integrating over Q, 


j sabe aed \ an : X(k,uo)duo 
(3.11) Go(k) I , tan Ri ea itu) 


This equation for the Fourier transform of the flux, p(r), is similar to equation 
(2.4), which represents the Fourier transform of the flux in the case of plane 
symmetry. However, there is one important difference, which is that, for the 
plane case, S(u) in equation (2.4) is obtained directly from the plane surface 
distributions, whereas X(k,xo) is obtained from an integral of the surface 
distribution, taken over the surface of the sphere. 

If we represent the integral on the right-hand side of equation (3.11) by 
Y(k) we obtain for the flux p(r) 

1 V(k)e “dV; 


(3.12) = (24)° J 1—(c/k) tan 'k 


where dV; is the three-dimensional volume element in k space. 
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Since Y(k) does not depend on angle and is an even function of k, we obtain 


"@ kY(k)sin(kr)dk 


(3.13) r= 3 A 1—(c/k) tank’ 


This integral is similar in form to that obtained from equation (2.4) and 
again the type of contour used in Fig. 1, with poles at k = +ixo and branch 
cuts from k = 7 to i~ and k = —i to —io, will be used in evaluation of 
equation (3.13). The contour will be taken either in the upper half plane or 
lower half plane as is necessary to ensure convergence of the integrals. The 
case ¢ <1 will be considered here and the case c = 1, when xo = 0, will be 
considered to be the results of taking the limit as c approaches unity. The case 
of ¢c > 1 is also considered in CHP. 

By the procedure which has led up to equation (3.13) we have replaced our 
two-medium problem by an infinite medium problem with a surface shell 
source at r = a. However, the solution to our problem still appears naturally 
in two parts; that is the solution for r < a and that for r > a. 

If we had the exact angular distribution on the surface of our sphere, we 
would have from equation (3.13) only a solution regular at the origin of the 
infinite medium source-free equation for r < a (CHP, Section 17). This solution 
would be of the form 


das(r) = A [sinh (Kor) /kor]. 


That is, the transient* part of our solution should vanish for r < a, and in order 
to satisfy the condition that the asymptotic solution should vanish for r<a, 
we must subtract off an amount of the solution of the infinite medium source- 
free equation to cancel ¢,,(r). A will be determined by 


(3.14) A= Pas (7) | pmo. 


For our approximate calculations we will not have the exact surface dis- 
tributions and consequently the solution evaluated for r < a, from equation 
(3.17) will have a small proportion of transient term mixed in. This will only 
be important for a sphere of small radius. However, this will mean that our 
solution for r <a will not vanish completely inside the sphere. We can, how- 
ever, subtract enough of the source-free solution to make our solution vanish 
at one point, the point r = 0 being the obvious choice. In other words we 
subtract off an amount of the source-free infinite medium solution, regular 
at the origin, given by A [sinh(kor)/xor] where A will be determined from 


(3.15) A = p(r) [rao 


where p(r) is given by equation (3.13). 
In practice we need not go through the contour integration just described, 
but we can make use of the results already given for a point source by CHP. 


*Strictly speaking in the case of a finite region (such as the interior of a sphere) the separation 
of p(r) into pas(r) and pnon-as(7) is not unambiguous. However, since (V?+x5)pas(r) = 0, the 
statement in the text should be understood in the sense that if we knew the source distribution 
exactly, then for r < a we would have (V?+x3)p(r) = 0. 
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This is done as follows. We rewrite equation (3.12) as 
a (c/k) tan~*k \ -as 
(3.16) p(r) = apf ees (c/k)tan"k e Y(k)dV,. 
Now the first term simply represents the direct unscattered contribution to 
the neutron flux from our surface source distribution. Also, 
(c/k) tan *k 
1- (c/k)tank 


is proportional to the Fourier transform of the isotropic point source solution 
(CHP, Section 14). This means that we can invoke the convolution theorem 
for Fourier transforms and write the flux in the following form 

(3.17) p(r) = Z(r)+ef pollr’—t|) Z(r') dV», 

where Z(r) is the direct unscattered contribution from the surface source and 
po(r) is the point source solution, 


a) oll elec 
; Po A4nr 0c 0 £ Mh ig . 


Z(r) is readily calculated for a specific surface distribution ¥(r’,Q)|,22, as 
follows, using Fig. 2, 


oe 


R 


Ae 


Fic. 2. The determination of the unscattered neutron density at r from a surface source 
onr =a. 


(3.19) F(P,Q) = e*Y(8',Q) | m0 


where f(r,Q) is the angular distribution at r due to the surface distribution 
at r = a. Then 


(3.20) Z(r) = s f f(r,Q)d2 = f f(r,Q)du, 


where pn, = Q.r/r. Thus 


(3.21) Z(r) -S searta 


where un, =Q.a/a and ¥(a,u,) is the surface distribution. 


Since 
a? = r?+ R?—2Rrp, 


and 
r? = a?+ R?+2Ran, 
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we can rewrite equation (3.21), taking R as the variable of integration, as 


1 r+a a ( — a8 
(3.22) Z(r) = oy ar v\ a, Ra 1- Re dR. 
a—r, if r<a 















Since we have obtained expansions in orthogonal polynomials for ¥(a,u,) for 
various values of the sphere radius, it suffices to calculate 
















rile rta eer -#) 
(3.23) Z,(r) = = Ee Ra 1——ar- dR 
a—T, itr<a 





corresponding to the various powers of yu” in ¥(a,u,). These components of 
Z(r) can then be combined according to the coefficients of u® in the expansion 
(3.7) to give the final approximation to Z(r). 

For example, we obtain 








\ 1) —Vr2-a? —(r+a 2 2,| Ee wy 
(3.24) Zr) = LU “ 1 ett) _ (72g | Babee) 
ete} 
(rta) JS’ ” ue 
a —(a—r)__ (rte) 72 | Bla) Meee 
7 at : =e) a-r r+a Af Sh 





The other Z,(r7) used in our calculations are listed in Appendix I. In the 
above formulae E,,(x) is one of the functions E,(x) defined by 






(3.25) E,(x) = f ou "du. 
1 





Some properties of these functions and tables for various m are given in CHP, 
Appendix A. 

Since po(r) as given by equation (3.18) is not known as a simple analytic 
expression, and since Z(r) is spherically symmetric, the integration in equation 
(3.17) involves a double numerical integration for each Z,(r). This can be 
reduced to a single integration if a suitable accurate expression for po(r) can be 
found. Such an expression is derived in Appendix II. 

All that remains now, for any particular problem is to evaluate (3.17), and 
then determine the amount of source-free solution to be subtracted, using 
(3.15). 

The asymptotic solution from (3.17) for r>a will have the form 













—gr 












(3.26) Pas. outside (7) = B r 






and subtracting off the amount of sinh(kor)/«or given by (3.15) the final 
asymptotic solution will be 






—Kgr 


(3.27) u(r) = Bi——A 





sinh(«or) 







Kor 
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and in the limit as c approaches unity, with xo approaching zero 






Pas(1) as 2A 





whence 





ae 
(3.28) A= Be 4: 








A is readily obtained from (3.15) and (3.17) as 


(3.29) A = Z(0)+42c¢ J oolr)Z(yrar. 
0 





B can be obtained from (3.17) by noting that in the integral only the term of 
the form 








e~colF—F'l /|p—r’| 





will contribute to B and this term can be expanded as 


e*\r—r'| _ sinh"**”? ¢~"*’ terms decreasing faster 
e*|r—r'| ee rengen toaile SS 
\r—r | kor 





(3.30) 


forr— ©, 








| than 





The terms decreasing faster than e~"°’/r do not contribute to the value of B, 
whence, since Z(r) is of order e~’ also, we obtain 







—cdKo sinh (xor’) 


(3.31) B= ac J “try 2" )rdr’, 






which can easily be expressed in terms of the Fourier transform of Z(r) (with 
xo replaced by —72k), which in turn is simply related to Y(&). 
We obtain finally 










—ca’ dks a. : e Od ny 
3.02) B= = ym) d Ba RT SRE) PED EEL ST 
(3.82) 2 dc J-1 wv (a,n)du -1 V[1+2koupo— xo (1—u*—u)] 





For c = 1, «xo = 0, 





(3.33) B= —3a" f pw(a,u)dp. 





Having now obtained the expression for A and B it remains only to consider 
the various values of a for which d is desired and to use the iterated angular 
distributions from (3.5) and (3.7) in the computation. 

We will consider the limiting cases of small a and large a first of all since it 
is only in these cases that accurate results are known for \ (Davison and 
Kushneriuk 1946; Davison 1951). These results are: 


.7104+0(1/a), 
4/3+0(a). 


For large a, the spherical surface of the sphere approaches a plane surface 















for large a, r 
(3.34) 






for small a, r 
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and we return to the Milne problem, which has already been solved by the 
iterative method proposed (Conkie 1959), the result being 


\ = .7113+0(1/a). 


For small a, we have the simple result that 


¥(a,u) = 494.) terms of order y 


and to this order is independent of angle. Thus 


(3.35) B = 322 +0@). 


Also for the interior of a small sphere, the probability of a collision is small, 
and consequently the first term in the expression (3.17) for p(r) dominates 
inside the sphere. This means, since it is p(0) which determines A, that only 
the Z(0) terms in equation (3.29) are important for small a, the scattered 
terms being of order a relative to the unscattered terms. In particular we need 


Z,(0) for small a. 
This is easily calculated from equation (3.24) from the formula for r<a, and 
gives 


Z.(0) = 2+0(a) 
and then 
(3.36) A = 2q0/a’. 
Substituting (3.35) and (3.36) into (3.28), we get, for small a 
\ = 4/3+0(a). 
Combining our iterative results for large a and small a, we have then 


for large a, \ = .7113+0(1/a), 


(3.37) 
for small a, A = 4/3+0(a). 


Comparing these results with (3.34), we see that for large a we have an 
error of about 0.1% and for small a we obtain the exact result. These results 
might lead us to hope that for intermediate a, high accuracy can be attained 
by this method. Unfortunately there are no very accurate results for inter- 
mediate a with which comparisons can be made. There do exist some approxi- 
mate results obtained from various methods and these will be discussed in the 
next section. 


4. COMPARISONS WITH THE RESULTS OF OTHER METHODS 


We will consider the results of two other workers for the black sphere 
problem, and we will use their results for comparison with the results of our 
iterative method. Among the results we will consider are those of Davison 
and Kushneriuk (1946), who used a combination of spherical harmonics 





308 : CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


calculations and perturbation calculations to estimate the graphical form of 
the dependence of A on a. We will also consider the variational results of 
Marshak (1947), who used a stationary expression for \. However, Marshak’s 
results were calculated using a rather crude trial function and we will not lay 
much stress on this comparison. Indeed, it would be desirable to repeat 
Marshak’s calculation, using a better trial function before going into any 
detailed comparisons with his variational results. 

The improved angular distribution was expanded in shifted Legendre 
polynomials as in equation (3.7) and terms up to and including Pf(u) were 
used in most of the calculations. Estimates of the effect of retaining higher 
order terms will be given later. The integration in (3.29) was done by Simpson’s 
rule with an interval of Ar’ = .02 for a range of values of a. (The integration 
was stopped at values of r’, where only further contributions of order 10-* 
were being made.) The integration was repeated with Ar’ = .035 for a = 1.0 
in order to assess the truncation errors due to the quadrature formula used. 
This comparison will be given below as well. 

The results obtained are given in Table III and the corresponding values 
from Davison and Kushneriuk’s estimated curve are given for comparison. 


TABLE III 
Values of \ as a function of sphere radius 





Davison’s Iterative 
value of X value of \ 
4/3 4/3 
1.256 1.2557 
1.195 1.1757 
1.145 1.1339 
1.100 1.0761 
1.060 1.0535 

. 988 . 9876 


Davison and Kushneriuk’s curve and a curve drawn from the iterative 
results are shown in Fig. 3. The curves differ by less than 2% from a = 0.4 
toa = 0.8, by about 1% at a = 1.0, and are very close at a = 1.5. The precise 
numerical value obtained from the estimated curve cannot of course be taken 
very seriously, but the results seem to show that the estimated curve follows 
quite well the accurate form of the curve for \ as a function of a, the greatest 
errors being in the neighborhood of a = 0.7-0.8. 

Marshak’s variational value is \ = 1.024 for a = 1.0, which is about 3% 
below the iterative value. Presumably this value would be improved con- 
siderably if a more accurate trial function were used. It should be noted that 
owing to the form of variational principle used, Marshak’s approximate value 
should indeed lie below the accurate value. 

Asa check on the internal consistency of the calculation, further calculations 
were carried out retaining only up to P}(u) and Pf(u) terms in the expansion 
(3.7) for a = 1.0. The results obtained for Z(0) for both of these distributions 


and for the Pf(u) distributions are given in Table IV. 
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--- DAVISON'S CURVE 
— ITERATIVE RESULT 
© MARSHAK'S RESULT 


RADIUS ,a 


Fic. 3. A comparison of the linear extrapolation length \, calculated by the iterative 
method, with other estimates of A. 


TABLE IV 


Values of Z(0) for various angular 
expansions, a = 1.0 


Order of expansion Z(0) 


. 72076 
. 72194 
72177 


We see that there is a very small difference between the Pj and Pf results, this 
difference being about .02%. Thus we can take the Pf results for Z(r) to be 
accurate enough for most purposes, assuming that the accuracy in Z(0) is an 
indication of the accuracy in Z(r). Any other quantity which we may wish is 
obtained from the second term on the right of equation (3.17), which involves 
an integration over the function Z(r), which should have an oscillatory error 
due to the orthogonal function expansion used. These other quantities should 
then be even more accurate than Z(r). 

For example, for the P}(u) distribution we obtain \ = 1.0514, which differs 
by only 0.2% from the Pf(u) value, A = 1.0535. 

By changing the interval in the Simpson’s rule integration from Ar = .02 
to Ar = .035, we obtain from the P}(u) distribution, \ = 1.0442, which differs 
by about 1% from the presumably more accurate value. 

Considering all these factors, and the fact that we have used an approxima- 
tion to po(r), our iterative results should be accurate to better than 1%, the 
errors lying mostly in the numerical approximations made in evaluating the 
various integrals involved. 
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5. EXTENSION TO GENERAL PROBLEMS WITH SPHERICAL SYMMETRY 


Having discussed the particular problem of the black sphere in Section 3 in 
some detail, we can easily determine how a general multilayer problem should 
be tackled. The iterative procedure follows through the first two steps and we 
obtain the iterated angular distribution at the surfaces r = r, from equation 
(3.4). It should be noted here that although in most problems the angular 
distribution will not vanish over part of its range, as it does vanish for 
Q.1r/r > 0 for the black sphere problem, it is still best to expand the distribu- 
tion in two sets of shifted Legendre polynomials, one adjusted to the range 
— 1<yu< 0 and the other to the range 0 < uw < 1. Then the coefficients in the 
angular expansion for u > Oare obtained by inserting the density and value of c 
for ri-1<7<™% in equation (3.4) et seg. and the coefficients for u < 0 are 
obtained by inserting the values of density and c for r; <r < r441 in equation 
(3.4) et seg. This takes into account the boundary conditions that at the 
surface r = r,, neutrons with » > 0 have come from r < r;and neutrons with 
u <0 have come from r > 7;. 

Now we go on to step 3 of our iterative procedure. For ry < r < ri41 only 
the surface sources at r = r; and r = 7,4; contribute to the density. The 
complete solution in this region is then obtained as a superposition of the 
contribution from each surface, and the amount of source-free solution which 
is added in is determined as before by p(r = 0) = 0 and p(r>ri4:) = 0, 
where p(r) is determined from the surface sources at r = r; and 7 = 7444. 

Going on in this way determining the contribution from the boundary 


surfaces of each region we can determine the iterated density for that region 
in the form of equation (3.17), whence one numerical integration will give us 
the explicit values for the density in the region being considered. 


6. CONCLUSION 

We see then, that our iterative method can readily be extended to problems 
with spherical symmetry. The calculations required are considerably more 
complicated than the original spherical harmonics calculations or even than 
the iterative calculations for problems with plane symmetry. For a preliminary 
design calculation on a spherical system the amount of computation and 
analysis required is probably prohibitive, and a P3; or Ps calculation will 
prove less troublesome. On the other hand for such a problem where the usual 
spherical harmonics method converges slowly or where more accurate results 
are required, the more elaborate scheme proposed here seems capable of giving 
results of extremely good accuracy. 
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APPENDIX I 
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APPENDIX II 
An ANALYTICAL EXPRESSION FOR THE POINT SOURCE SOLUTION 


We wish to obtain a simple, accurate expression for 


1 —rim Oi 
(11.1) pir) = po -f elem" Ss 


where 
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(CHP, Section 14). Since g(c,u) is a function of uv? only, we will try and find 
an expansion of g(c,u) in the orthogonal polynomials p, (u?), such that 
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where 6nm is the Kronecker delta. 
The three term recurrence relationship between these polynomials is easily 


determined to be 
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Expanding 
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we obtain 
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The moments (u*") of g(c,u) are given in CHP (p. 78 et seg.) for various 


values of c. 
Let us consider c = 1.0 and N = 2 in the expansion (11.4). Then we obtain 












g(cu) = .8—.6859 p,(u?) —.4484 po(u?) 







.9902 — .3016 uw? —.4484 yt. 





Inserting this into equation (II.1), we obtain 
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where the £,(x) are the same functions referred to in Section 3. 
Comparing this with the tabulated values given by CHP, we obtain the 
results given in Table V for e(r) = 4are"p,(r). 







TABLE V 


Comparison of N = 2 approximation with exact results. 
c=1.0 

















n 


€ exact € approximate % error 













0 1.0000 . 9902 -1 

5 . 8460 . 8479 +.2 
1.0 .7612 . 7635 +.3 
2.0 .6568 .6571 Very small 









We see that the error in €(7) is greatest at r = 0, and drops off fairly quickly, 
being oscillatory in sign. This latter is a useful property since it will mean 
some cancellation of errors in integrals involving our approximate function. 
The agreement between the two sets of results is good considering the rather 
low order of approximation. 

It was felt that the NV = 3 approximation would be good enough for the 
calculations done on the black sphere, and the results for this case are as follows. 










g(c,) = 1.0055 —.6220 u2+.5130 u!—.7050 us 






and 
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(11.7) pi(r) = 24 1.0055 .62208.(r)+.5130E(r)—.7050E4(r)t. 






Table VI shows the comparison between the e(r). 
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TABLE VI 


Comparison of N = 3 approximation with exact results, 
c=1.0 


€ exact € approximate % error 


1.0000 1.0055 +.5 
. 8460 . 8450 -.1 
. 7612 . 7604 -.1 
. 6568 .6570 Very small 


The results show that this approximation should be accurate enough for 
most practical purposes although it should be noted that it is very easy to 
increase the order of approximation in the expansion of g(c,u). 

Although not used in our calculations, the results of such an expansion for 
c = 0.5 may also be of interest, and will indicate the usefulness of this expan- 
sion for other values of c. 

We obtain, for the V = 2 approximation, 


(11.8) sitet ZL} 66+ s22ans(e)-+.20188,()}. 


Table VII shows the comparison with exact results. 


TABLE VII 


Comparison of N = 2 approximation with exact results. 
c=0.5 








€ exact € approximate % error 


1.0000 . 9866 ie 
1.1529 1.1551 i 
1.2412 1.2454 
1.3504 1.3521 





Again the error is oscillatory and largest near r = 0. A low order approxima- 
tion should be adequate for calculations with this value of c, and, the indica- 
tions are, for all values of ¢c of physical interest. 





ON THE PRODUCTION OF Be’, Mg’, AND Ni*®* IN 
THE SLOW NEUTRON FISSION OF U**! 


J. C. Roy 


ABSTRACT 


The yields of 53-day Be’, 21.3-hour Mg*, and 56.6-hour Ni® from the slow 
neutron fission of U?* have been investigated. Upper limits of 3X10-7% and 
4.2X10-°% have been set for the fission yields of Be? and Mg** respectively. A 
fission yield of 2.0+1.0X10-8% has been found for the formation of Ni®. These 
results are compared with the current knowledge of the frequency of triple 
fission of U5, 


INTRODUCTION 


The possibility of slow neutron fission of U** in three fragments has been 
studied by many investigators with a variety of techniques. From the findings 
of these studies, reviewed by Whitehouse (1952), Perfilov (1958), Demers 
(1958), and Hyde (1960), some features of the ternary fission of U?* by slow 
neutrons have been well established while others are not yet clarified. 

One type of ternary fission, namely the emission of long range alpha particles 
in coincidence with two heavy fragments has been completely proved. The 
frequency of this process for slow neutron fission of U** is one in 325+100 
binary events. 

The emission of short range light particles associated with two heavy frag- 
ments is often described as another type of ternary fission. Cassels et al. 
(1947), Allen and Dewan (1951) working with counters, Green and Livesey 
(1947, 1948), and Titterton (195la), working with nuclear emulsions have 
found that one slow neutron fission of U?* in about 100 resulted in the emission 
of short range particles with a charge of 1, 2, or slightly higher. Tsien and his 
collaborators (1947, 19475) also observed short range particles of about the 
same abundance in nuclear emulsions, but they could not decide on the 
source. They may result from fission or they may be protons and other light 
particles recoiling from collisions between the nuclei in the emulsion and 
fission fragments. Marshall (1949), Demers (1953, 1958), and Perfilov (1958) 
support the latter view. deLaboulaye et al. (1953, 1954) using cloud chambers 
concluded that if the emission of short range light particles exists at all, its 
frequency is smaller than 1+3 in 1000 binary events. 

The possibility that some of these fragments are beryllium nuclei have been 
investigated by Cook (1952), who set an upper limit of 10-°% for the yield 
of Be’ and by Flynn et al. (1956), who set an upper limit of 4X 10-*% for the 
yield of Be!®. Titterton (195la, 6) from the examination of 250,000 fission 
tracks has observed no instance of Li’, B’, or Be’. Recently Albenesius (1959) 
has reported evidence for the formation of tritium in the ratio of one to 10¢ 


1Manuscript received October 28, 1960. 
Contribution from the Research Chemistry Branch, Atomic Energy of Canada Limited, 


Chalk River, Ontario. 
Issued as A.E.C.L. No. 1164. 


Can. J. Phys. Vol. 39 (1961) 





316 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


binary fissions. He interpreted the formation of tritium as due to ternary 
fission. From these studies one may conclude that the frequency of emission 
of a short range particle is not yet established, with the exception, perhaps, 
of tritium. 

In the present work a further radiochemical search for Be? was made. 
Radiochemical investigations cannot show what type of fission produces Be? 
and therefore the discussion of possible mechanism will be postponed until 
the formation of Mg?® and Ni® is discussed. 

A third type of ternary fission, the splitting of U®** in three fragments of 
comparable mass, can be differentiated and is usually called tripartition. 
Theories of nuclear fission (Bohr and Wheeler 1939; Present and Knipp 1940; 
Present 1941) do not rule out tripartition but say little about the frequency 
of the process. 

Tsien and his collaborators (1947, 19476), Chastel and Vigneron (1949), 
Perfilov (1958), and Catala et al. (1959, 1960) have observed triple tracks in 
nuclear emulsions which they ascribed to tripartition. Using the same tech- 
nique Demers (1946, 1947, 1953, and 1958), Green and Livesey (1947), Titter- 
ton (1951a), and Dutta (1953) upon examination of 1500, 25,000, 250,000, 
and 20,000 fission tracks respectively could not definitely assign any of the 
observed events to tripartition of U***, Rosen and Hudson (1950), from a 
careful study of tripartition of U*** with triple ionization chambers, concluded 
that triple fissions occur to the extent of 6.7+3.0 per 10° binary fissions. 

In the wartime radiochemical research on fission, attempts were made to 
detect radioactive products in the mass range of 30 to 60 which could have 
been formed by tripartition, but the results were negative (Metcalf e¢ al. 
1951). Upper limits of about 10-*— 10-°% were set for the formation of such 
products. The findings of these studies on tripartition are collected in Table I. 

The frequency of tripartition obtained by Perfilov and Tsien et al. seems to 
be too high in comparison with that obtained by the other investigators. The 
first group of workers may have mistaken for a tripartition a triple track 
due to a binary fission event plus a heavy recoil originating in the emulsion at 
approximately the point of fission (cf. Demers 1958 for a discussion of the 
question). Actually Dutta (1953) observed one tripartition, but he ascribed it 
to the interaction of a fast neutron of cosmic origin. Demers (1958) does not 
exclude the possibility that the triple events observed by Rosen and Hudson 
(1958) could be due to collisions of binary fragments with C, N, or O atoms 
at the point of fission. 

In summary, it appears from the work of Titterton, Catala e¢ a/., and Rosen 
and Hudson that the frequency of tripartition is one or less per 210° binary 
fissions. Physical methods are able to detect single events whereas the yield 
of radioactive product represents a fraction only of the total tripartitions. If 
one assumes that radioactive products in Table I form about 1% of the total 
tripartition, an upper limit of about one per 10‘ binary events is obtained 
from the radiochemical studies. 

Another type of fission which could produce the radioactive products 
investigated in our study is the quadripartition of U* by thermal neutrons. 
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TABLE I 
Summary of the findings on the tripartition of U** by thermal neutrons 


Particle masses 


triple fission 


Authors M, M:2 M; Ratio = binary fission 
(A) Photographic emulsions 
Perfilov (1958) 62 113 60 1/5000 
73 76 86 
Tsien et al. (1947, 19476) 127 te 32 1/6000 
2 
Catala et al. (1959, 1960) be a " 1/175,000 
Demers (1958) None in 1,500 b.f. <1/1500 
Green and Livesey (1947) None in 25,000 b.f. <1/25,000 
Titterton (1951a) None in 250,000 b.f. <1/250,000 
Dutta (1953) None in 20,000 b.f. <1/20,000 
(B) Ionization chambers 
Rosen and Hudson (1950) 6.7/10° 
(C) Radiochemical studies 
Metcalf et al. (1951) 
Nuclide Yield 
87-day S® <1.8X10“% 
160-day Ca* <1.4X107-'% 
4.7-day Ca“ <10°% 
Scandium isotopes <10-'*% 
45-day Fe? <4X10°% 








Tsien and his collaborators (1947, 1947a, 6, c; Ho et al. 1946; Chastel and 
Vigneron 1949) claimed to have observed cases of quadripartition with a 
frequency of one in 5000 binary fissions. However, Perfilov, who has observed 
similar cases, attributed them to elastic collisions of each fragment with nuclei 
of elements in the emulsion. Demers (1958) supports the latter view. This is 
also in agreement with the observations of Green and Livesey (1947), who 
found no quadripartition in 25,000 events, and of Titterton (1952), who saw 
only two tracks in 600,000 events which could be possibly classed as quadri- 
partitions. 

One may conclude from this review that additional evidence is required to 
establish with certainty the occurrence of tripartition and quadripartition 
and that the frequency of such events is low, in particular the yield of any 
product resulting from these processes is probably less than 10-°%. For such 
a low yield a search for suitable radioactive products seems to be the best 
approach to explore anew tripartition or quadripartition of U* by slow 
neutrons. Identification of a radionuclide which cannot be produced by 
reactions other than fission should provide evidence for the processes. Of 
course the formation of such a product by binary fission should also be ruled 
out. Three radionuclides, 53-day Be’, 21.3-hour Mg?*, and 56.6-hour Ni® are 
the best candidates for such an investigation. These three nuclides cover a 
wide range of mass; in a thermal neutron irradiation of uranium they can 
be formed only by fission, except for very rare processes like double neutron 
capture; they have suitable radioactive properties and they can be easily 
purified. The present paper reports the results of a search for these three 
radionuclides in the slow neutron irradiation of U**. 
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EXPERIMENTAL 
(a) Irradiation and Purification 

One sample of natural uranium trioxide (~150 mg) and one sample of 
reactor grade uranium dioxide (~5 mg) enriched to 938% in U?* were used 
as starting materials. These were sealed in a quartz tube and irradiated in an 
iron capsule for 16.6 hours in an empty fuel rod position in the NRX reactor 
(Hurst 1955). The reactor neutron flux in the position used was about 
7X10" n/cm? sec; about 94% of the flux is thermal, 3% isin the resonance 
neutron region, and the other 3% are fast neutrons (Roy and Wuschke 1959). 
After the end of the irradiation, the iron capsule was taken out of the reactor 
by means of a pneumatic carrier and received and opened in a shielded cell 
where simple remote controlled operations could be done. The contents of the 
capsule were poured into a beaker containing known amounts of Be, Mg, Ni, 
and Ba carriers. The barium was extracted in order to determine the number 
of fissions. The irradiated uranium oxide was heated for about half an hour 
in concentrated HCl; this treatment dissolved about one-half of the oxide. 
The solution was transferred to another beaker, diluted with water, and the 
barium was precipitated as the sulphate. The BaSO, was filtered and then 
put aside for further purification. The filtrate was evaporated to a small 
volume, made 10 to 12 N in HCl, and passed through a Dowex-1 ion exchange 
column; this step eliminated the uranium and some fission products (Kraus 
and Nelson 1955). The eluate was evaporated to a small volume, neutralized 
with ammonium hydroxide, and beryllium hydroxide was precipitated by the 
addition of a few drops in excess. The hydroxide was filtered and kept for 
further purification. Sodium hydroxide was added to the filtrate and, on 
boiling, nickel and magnesium hydroxides precipitated. The precipitates were 
filtered. The hydroxides were dissolved in HCI and the solution was passed 
twice through Dowex-1 in 10 N and 1 N HCI respectively. Ferric hydroxide 
was precipitated in the solution several times and the nickel was separated 
from the magnesium by a dimethylglyoxime precipitation followed by a 
centrifugation. Complete purification of the nickel was achieved by the 
procedure described in Kleinberg (1958). Finally nickel was electroplated on 
a weighed platinum foil, weighed, placed on an aluminum tray, and its 
radioactivity determined with an anthracene beta counter described in the 
next section. 

The magnesium fraction in the supernatant resulting from the dimethyl- 
glyoxime precipitation was boiled with sodium hydroxide. The magnesium 
hydroxide precipitate was dissolved in concentrated HNO; and two MnO, 
scavengings were done. The magnesium was precipitated as the ammonium 
phosphate by the method of Epperson (Treadwell and Hall 1945). The precipi- 
tate was dissolved in 1.5 N HCl and passed twice through Dowex-50 at 
a temperature of 50° C (Milton and Grummitt 1957). Finally the magnesium 
was precipitated as the ammonium phosphate, filtered on a ‘‘Millipore’’ filter 
paper, dried at 50° C, weighed, placed on an aluminum tray, and counted 
with the anthracene beta counter. 

The purification of beryllium hydroxide obtained in the first stage of the 
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separation was done according to the method of Buchanan (1958). Finally 
the beryllium was transformed to BeO, deposited on an aluminum tray, 
weighed, and counted with the Nal(TI) scintillation spectrometer described 
in the next section. 

The BaSQO, precipitate obtained at the beginning of the separation was 
dissolved in a sodium carbonate solution. The purification was done according 
to the procedure described in Kleinberg (1958). The final BaCrO, precipitate 
was dissolved in dilute HCI and the solution made up to a known volume. 
Aliquots were pipetted out on thin plastic films and the activities of the 
sources were measured by 44% 8—y coincidence techniques (Campion 1959). 












(b) Counting 

The sources of Mg(NH,)PO,.6H:O and Ni were covered with 42 mg/cm? 
of aluminum. This eliminated several counts of long-lived unidentified activi- 
ties: it also completely cut out the low energy beta particles of Ni®* and most 
of those from Mg?’. Thus the anthracene counter detected only the beta 
radiations of Al?* and Cu® in equilibrium with their parents Mg?* and Ni®, 
The distance between the source and the detector was about 1 cm of air. The 
counting efficiency of the apparatus for the beta radiations of Al*?* and Cu® 
was 28%. The background of the detector was 8.5+0.2 counts/min. 

The intensity of the 0.48-Mev gamma ray of Be’ in the BeO source was 
measured with a 3in.X3in. Nal(Tl) crystal used in conjunction with a 
100-channel pulse height analyzer. The efficiency of detection of the crystal 
for the 0.48 gamma ray was 11%. 
















RESULTS 
(a) Determination of Fission Yields 

The fission yield of Be, Mg, and Ni were calculated from the saturation 
activities of the species and of Ba'*® produced in the same sample using the 
equation (Coryell and Sugarman 1951) 


(1) 


where Yj, is the yield of the species to be calculated, 
Sp is the saturation activity of the species at end of irradiation, 
S, is the saturation activity of Ba'° at end of irradiation, 

and_—“‘Y,_ is the yield of Ba"® taken to be 6.44% (Katcoff 1958). 



















The saturation activities of Ba'° were obtained from the measured equilibrium 
activities of Ba'°-La'®, with the application of a correction factor (2.14) for 


transient equilibrium. 


(b) Beryllium-?7 

In the examination of the two sources of BeO with the scintillation spec- 
trometer no evidence was found for the presence of the 0.48-Mev gamma ray. 
Thus upper limits only could be set for the formation of Be’; it was done in 
the following manner. Alternate measurements of the activities of the source 
and of background were repeated several times for periods of 30 minutes (for 
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sample 1) and of 10 minutes (for sample 2) over the energy range of 0.05 to 
2 Mev. The absolute difference of the activities between the source and back- 
ground in the 0.42 to 0.54 Mev interval was taken as a measure of the limit 
of the detection of the 0.48-Mev gamma ray of Be’. The mean of the values 
of the several measurements recorded in succession was taken as the best 
estimate for the limit of detection of Be’. The data for the calculation of the 
fission yield of Be? are given in Table II. In the calculation of the yield, the 


TABLE II 
Fission yield data for 53-day Be’ 


Run 1 Run 2 


Nature of the uranium target Natural Enriched 
Maximum activity of Be’ extrapolated to the end of the 

irradiation (counts/min) t2 0.9 
Abundance of the 0.48 Mev y-ray (%) 11 11 
Chemical yield (%) 73.8 ‘ 
Counting efficiency (%) 11.0 
Disintegration rate of Be’ at the end of the irradiation 

(d.p.s.) 2.23 
Correction of disintegration rate for destruction of Be? 2.45 
Fractional saturation (16.6 hours) 9.041073 
Saturation activity of Ba’ (d.p.s.) 1.85X10° 
Saturation activity of Be’ (d.p.s.) 2.71 X10? 
Fission yield (%) <9.4 X107 





destruction cross section of Be’? by thermal neutrons, the cross section of 
which is 4.810! barns (Hanna 1955), has been taken into account. 


(c) Magnesium-28 

The source of Mg(NH,)PO,4.6H,0 had a counting rate of 10 counts/min 
above background 31 hours after the end of the irradiation. This activity had 
decayed to 8.4 counts/min 4 days after the beginning of the counting of the 
sample. Since it represented 4.5 half-lives of Mg?’, 1.6 counts/min was taken 
as the maximum possible amount of activity due to Mg?’ 31 hours after the 
end of the irradiation. Examination of the source with the scintillation spec- 
trometer gave similar limits for the gamma radiations of the doublet Mg?*—Al?®. 
The value of 4.210-°% is considered an upper limit for the fission yield of 
Mg?’ since no positive identification of Mg?® was made. The data for the 
determination of the fission yield of Mg?’ are presented in Table III. 


TABLE III 
Fission yield data for 21.3-hour Mg?8 





Nature of the uranium target Natural 
Maximum activity of Al? daughter at the end of irradiation (counts/min) 4.0 
Chemical yield (%) 41.3 
Counting efficiency of the 2.87-Mev beta-ray (%) 28.0 
Disintegration rate of Mg** at the end of irradiation (d.p.s.) 0.58 
Fractional saturation (16.6 hours) 0.463 
Saturation activity of Ba‘? (d.p.s.) 1.85109 
Saturation activity of Mg*® (d.p.s.) 1.2 
Fission yield (%) <4.2X107-9 
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(d) Nickel-66 

Evidence for the formation of Ni®* in the slow neutron fission of U?* has 
been found. An activity which decayed with a half-life close to that accepted 
for Ni®® (56-hour from the compilation of Strominger et a/. 1958) has been 
measured in the two sources of Ni. Decay curves of the total activity and 
of that attributed to Ni® are drawn in Fig. 1. From sample 2, the more 
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Fic. 1. Activities of Ni® fraction: curves 1 and 3 gross decay curves; curves 2 and 4 
activities attributed to Ni®. Residual activities of 16 and 1.9 counts/min have been subtracted 
from curves 1 and 3, respectively, to obtain 2 and 4. 


active source (from enriched sample), it was shown that the activity was 
definitively Ni®* by separation and measurement of the 5.1-minute daughter 
Cu®*, Four successive separations of copper sulphide in 2 N HCl from the 
purified nickel solution gave a CuS fraction with an average activity of 
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12 counts/min, 5 minutes (or one half-life) after the end of the separation. 
The activity of Ni®® from the decay curve 2 (Fig. 1) is 26 counts/min at the 
time of separation. This established clearly that the measured activity (curves 
2 and 4, Fig. 1) was Ni®. 

The only reactions other than neutron fission of uranium which could 
possibly produce Ni®* are successive neutron capture in Ni®™ as follows: 


Ni*——(n,y)——> Ni®—— (u,v) — > Ni 
stable 2.36-hour 56.6-hour 


(1.0% of natural Ni; cross sections 2 barns, 20 barns (Hughes and Schwartz 
1958)). 

The amount of Ni required to form the observed activity is calculated to be 
0.1% in the natural uranium sample and 10% in the enriched sample. Analyses 
of the two samples gave 5 and 50 p.p.m. of Ni in natural and enriched samples 
respectively. Thus the reactions given above cannot be responsible for the 
formation of Ni®®. The data on the fission yield of Ni®® are assembled in 
Table IV. Considering the very low counting rate, which introduces a large 


TABLE IV 
Fission yield data for 56.6-hour Ni® 


Run 1 


Nature of the uranium target Natural 

Activity of 5.1-minute Cu® daughter extrapolated to end 
of irradiation (counts/min) 

Chemical yield (%) 

Counting efficiency of the 2.6 Mev beta of Cu® (%) 

Disintegration rate of Ni® at the end of irradiation 

Fractional saturation 

Saturation activity of Ni® (d.p.s.) 

Saturation activity of Ba" (d.p.s.) 1.85X 10° 

Fission yield (%) 2.9 <X10-° 


won 
OS™ OO StEt 
cco 


Neo 
an 


error in the measurement, the agreement between the two yields is reason- 
able. More weight is given to experiment 2 and the best value for the fission 
yield of Ni®® is taken to be 2+1X10-8%. 

The nature of the process which resulted in the production of Ni® is con- 
sidered in the discussion which follows. 


DISCUSSION AND CONCLUSIONS 


A summary of the results are given in Table V. Unfortunately, because of 
time requirement, it was not possible to dissolve all the uranium target. It 
could then be argued that some of the elements could have been preferentially 
leached from the target. If this is the case, a systematic error amounting, 
perhaps, to a factor of 2 would arise in the yields given in Table V. 

The present result for Be? reduces, by a factor of 30, the limit of 10% 
set by Cook (1952). It shows, in agreement with other investigators (Flynn 
et al. 1956; Titterton 195la, b) that the formation of beryllium nuclei in 
the thermal neutron fission of U?* is an extremely rare event. 
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TABLE V 
Summary of the results 






Yields (%) in the slow neutron fission of U% 


Sample Be? Mg?8 Ni® 


1 <9.4X107 <4.2X10-9 2.9X10-° 
2 <3X107 — 1.9X10-8 
Best value <3X107 <4.2X10-* = 2.0+1.0X10-8 


















The upper limit of 4.2X10-°% for the yield of Mg?* indicates that the 
formation of fission products in the region of mass 30 are so rare that their 
detection will require major revision and refinements of current techniques. 
This result does not support the observation of Tsien et al. (1947, 1947); 
cf. Table |), who have reported one case of tripartition with a fragment of 
mass 32 in nuclear emulsions. 

It can be shown that fast fission of U* or U8 cannot account for the 
observed yield of Ni®*. The fast fission cross section of U?* and U*38 are 1.3 
and 0.31 barns respectively and the fast neutron flux in the irradiation position 
used in the experiment is about 3% of the thermal flux. Hence the fast fission 
yield of Ni®® in the enriched and natural sample would have to be 14,000 
and 300 times higher than its thermal fission yield in order to compete with 
thermal fission of U***, The yields of fission products (Katcoff 1958) indicate 
that fast fission in a reactor could possibly increase the yield of an element 
like Ni®® only by a factor of 3 to 5. It is therefore concluded that Ni®® is a 
product of the thermal neutron fission of U**, but the exact process cannot 
be decided on the basis of the radiochemical evidence. Binary fission cannot 
be ruled out because an extrapolation of the fission yields of light mass frag- 
ments (Katcoff 1958) gives 10-8-10-'""% at mass 66, and this range borders 
on the observed value. One can only conclude that the yield of Ni®* is equal 
to or smaller than 2.0+1.0X10-°% either for binary, triple, and quaternary 
fission. 

Assuming that Ni® is produced by tripartition it is of interest to compare 
the frequency with that observed in other investigations of tripartition. To 
do so it is necessary to further assume that Ni® is formed in 1% of all tri- 
partitions. This assumption is reasonable if Ni represents most of the yield 
of the mass 66 chain, if the splitting of the uranium nucleus in three frag- 
ments of comparable mass is the most favorable mode of tripartition, and 
finally if tripartition is not too narrow. On this basis the frequency of tri- 
partition is one in 5X10’ binary fissions, a value much less than those reported 
by other workers and shown in Table I. The discrepancy could be explained 
if the yield of Ni® in tripartition was very low; it would have to be 5X10°% 
to agree with the lowest value reported previously. Some causes which could 
have led to a confusion of a triple event for a tripartition in nuclear emulsions 
and in the work of Rosen and Hudson (1950) have been discussed in the 
Introduction. Our results indicate that the problem requires new measure- 


ments. 
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The nuclear emulsion techniques offer little possibility; it would take 
Herculean labors to examine enough fission events to settle the question. 
However, the type of experiment of Rosen and Hudson could be repeated 
with special attention to eliminate the interference of C, N, and O atoms 
(Demers 1958) in addition to all the other precautions taken by the same 
authors. 
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THE STATISTICAL BEHAVIOR OF INSTABILITIES IN 
DISPLACEMENT PROCESSES IN POROUS MEDIA! 


ADRIAN E. SCHEIDEGGER? AND EDWARD F. JOHNSON?® 


ABSTRACT 


The problem of growth of instabilities (fingers) in displacement processes in 
porous media is analyzed from a statistical viewpoint. The’ relative area 
occupied by fingers is represented as a ‘‘saturation”’, and the equations of motion 
for this ‘‘saturation”’ are derived. It is shown that these equations are analogous 
to the Buckley—Leverett equations of immiscible displacement, with a fictitious 
“relative permeability” being introduced. The latter can be calculated and thus 
the statistical equations of motion of a fingered-out front can be written down 
explicitly. These equations of motion can then be solved by the well-known 
method of characteristics. It is shown that the statistical theory does not lead 
to any stabilization of the fingers. 


1. INTRODUCTION 


The problem of instability has been of great importance with regard to 
displacement experiments in porous media and displacement processes in 
underground rock strata in which the displacing fluid is less viscous than the 
displaced fluid. It has been observed that, under the conditions just men- 
tioned, protuberances may occur which may shoot through the porous medium 
at comparatively great speed. This phenomenon is called “fingering’’. A 
typical example of a fingered-out displacement front (after Terry et a/. 1958), 


is shown in Fig. 1. 


0.15 pore 
volumes injected 


Fic. 1. Example of a fingered-out front after Terry et al. (1958). 


In previous theoretical investigations (Chuoke ef al. 1960; Scheidegger 
1960a, 6, 1961; Perrine 1960), the problem of fingering has been treated from 
a ‘“‘microscopic’”’ standpoint. In this, the displacement front was assumed as 
uniform, small disturbances were superimposed, the differential equations 
governing the growth of the latter were linearized, and the further behavior 
of the front was obtained by an integration. Any theory of this type is valid 
only for infinitesimally small (‘‘microscopic’’) fingers. The result is invariably 

1Manuscript received September 29, 1960. 

Contribution from Imperial Oil Limited, Calgary, Alberta, and Jersey Production Research 
Company, Tulsa, Oklahoma (U.S.A.). 

2Imperial Oil Limited, Calgary, Alberta. ; 

3Jersey Production Research Company, Tulsa, Oklahoma. 
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that fingers, once started, should grow indefinitely (this is simply a conse- 
quence of the linearization of the equations). Emphasis has therefore been 
directed towards a determination as to whether instabilities of a certain 
wavelength do become started at all. 

However, one of the big questions is whether an array of fingers could 
stabilize itself while it is progressing through the porous medium. In order 
to investigate this problem, it is not sufficient to deal with the linearized 
displacement equations since any existing fingers can no longer be regarded 
as infinitesimally small. Therefore, a different approach must be sought after. 

In the present paper, it is proposed to treat the fingers statistically. Thus, 
only the average cross-sectional area occupied by the fingers is taken into 
account, the shape and size of the individual fingers are disregarded. The 
displacement of one fluid by the other is regarded as complete with no mixing 
occurring at the finger tips. The relative cross-sectional area occupied by the 
fingers is then analogous to a fictitious “‘saturation’’ of the displacing fluid 
in the porous medium. It turns out that the equation of motion for this 
“saturation”’ is similar to the Buckley—Leverett equation for the saturation 
change in immiscible, stable displacement, with a special expression for the 
relative permeability which becomes, in the present theory, a “‘fictitious”’ 
relative permeability. It is then necessary to introduce a theoretical model 
for the pressure distribution during the displacement. In particular, it is 
assumed that the pressures in adjacent elements in the two fluids are the 
same. This simplification will result in a description of a limiting case. With 
this theoretical model, it is possible to obtain a theoretical expression for the 
fictitious relative permeability. The solution for the spreading of fingers in 
a linear displacement experiment with fingering is then obtained and it is 
shown that the present theory does not yield any possibility of stabilization 
of the fingers. 

Owing to the simplicity of the theoretical model employed, more than a 
qualitative agreement with experimental data ought not to be expected. In 
any theoretical treatment of fingers, simplifications either at the microscopic 
(linear wave theory) or at the macroscopic level (present theory) have to be 
introduced. It is believed that the present attempt is the first instance of a 
treatment of fingers on the macroscopic level. 


2. THE MUSKAT MODEL 


The present treatment of the fingering process will be based upon what 
has been called the Muskat model of displacement theory, indicating its 
originator. This model assumes that there is a sharp front between the dis- 
placing and the displaced fluid, and that on each side the flow proceeds 
according to Darcy’s law, based on the corresponding viscosity. The Muskat 
model does not specify whether the fluids are miscible or immiscible, since 
only saturations (concentrations) of zero or unity are allowed. 

The Muskat model represents an oversimplification of the actual facts 
inasmuch as some mixing will always occur at the interface between two 
fluids. Nevertheless, it is known that it yields acceptable results in many 


cases. 
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3. STATISTICS OF FINGERS 
As noted in the Introduction, in endeavoring to describe the fingering 
process from a macroscopic standpoint, one is concerned only with the average 
behavior of the two fluids involved. 
Thus, consider a linear displacement experiment parallel, say, to the 
+x-direction which progresses with time (#). For certain values of x and ¢, 


there will be some fingers present. 
The fingers present in the porous medium can be schematically drawn as 
represented in Fig. 2. The problem, then, is to set up a determining equation 


DISPLACEMENT > 


LEVEL" x" 


Fic. 2. Schematic representation of fingers. 


for the function 
(3.1) Si = 54(x,t) 


where s; represents the average relative area occupied by the 7th fluid at the 
level x. 

This is done as follows: Using Darcy’s law for each of the phases, the 
pore velocities v; and v2 (referring to the displacing and displaced fluids, 
respectively) are obtained according to the following equation of motion 


k 
(3.2) Se grad p, 


where k is the permeability of the medium, P its porosity, and yw, we, are the 
viscosity of the displacing and displaced fluid, respectively. Furthermore, p 
is the local pressure. 

The theory cannot be carried further at this stage without assuming that, 
at the level x, the pressure » does not vary along any line on that level (a 
modification could be made, however, allowing for a capillary pressure differ- 
ential between the two fluids). Thus at the level x, the pressure is uniform. This 
leads to an expression for the flow of the phase 7 across the unit of area 


(3.3) a= Sts grad p. 
Mi 
The problem is determined, if in addition to the equations of motion (3.3) 


the continuity conditions are considered 


a 
(3.4) at = div qi, 
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(3.5) Sits, = 1. 


It is assumed, then, that the fluids involved are incompressible. 
Equations (3.3), (3.4), and (3.5) constitute the equations describing the 
average behavior of fingers in the present theory. 


4, FICTITIOUS RELATIVE PERMEABILITY 


The statistical treatment of the fingers outlined in the last section leads to 
a description of the displacement process which is formally identical with the 
description of immiscible displacement by means of ‘‘relative permeabilities’’ 
(denoted by k;,). 

One could have arrived at such a description of fingers also from a direct 
analogy with relative permeability. Accordingly, by again considering a linear 
displacement parallel to the x-direction and an analogy with immiscible dis- 
placement, it stands to reason that the behavior of the “‘saturation”’ s (i.e. 
the relative area occupied by fingers) as a function of x and ¢ can be treated 
in a manner analogous to the saturation in the Buckley—-Leverett theory of 
immiscible displacement. 

Thus, a ‘fictitious relative permeability” is introduced which describes the 
seepage flow q; of each phase 


(4.1) qa = ~h tenes 
Mi 


where k is the total permeability, u; the viscosity of the phase under con- 
sideration, and grad p the pressure gradient. This, together with the con- 
tinuity equations 

as, 


(4.2) = —div qi, 


(4.3) Zz se 1, 

(where P is the porosity) fully describes the flow process if the two fluids are 

taken as incompressible. It is assumed that the &,’s are functions of s; only. 
It is well known that the above equations have an eliminant 


Os1 
ot 


with g, denoting the filtration velocity in the x-direction and 


P Os 
++? (s1)ge 5. = 0 


(4.4) P 


x bs Ri/u1 
aR} a ki/urtRe/pe2 


The problem now is to find a suitable expression for k;(s). From the analogy 
with “relative permeability’”’ only, it is not possible to arrive at such an 
expression. However, by recalling the remarks on the statistical description 
of fingers made in the previous section, it is apparent that one must set 


(4.6) ki = Si 
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Fic. 3. Fictitious relative permeability curves. 


The corresponding ‘‘relative permeability curves’ are shown in Fig. 3. Using 
eq. (4.5) and eq. (4.6) the expression for r as a function of s becomes 





<! a Si/p_ ae me 
(4.7) r(si) = S1/mi+(1—s1)/u2  m—141/s 


with 
(4.8) m = uo/m. 


It will now be convenient to drop the index 1 referring to the displacing 
fluid. Thus r and its first and second derivatives with respect to s become 
m 


(4.9) r(s) = tetra 


i al 
(4.10) r'(s) [s(m—1) 41)" 
2m 
[s(m—1)+1] 
The shapes of the r(s) and r’(s) curves for a typical case are those given in 
Figs. 4 and 5 respectively. It is obvious that r(s) has no inflection point, 
r’(s) no extremal values, so that no shock front can develop (cf. the discussion 
in Scheidegger’s (1960c) book, p. 224). 


(4.11) 7'(s) = 3(m—1). 
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5. THE SPREADING OF FINGERS 


In order to determine the behavior of fingers as described in the present 
theory, it will be necessary to find appropriate solutions of the eliminant 
equation (4.4). Solutions of this equation can be obtained most easily by 
the method of characteristics. The characteristic equations corresponding to 


(4.4) are 
6.1) o = 1'(s)ge(t)/P 


(5.2) ds/dt = 0. 


Thus, eqs. (5.1) and (5.2) describe the paths along which a constant saturation 
moves. Inserting eq. (4.10) into eq. (5.1), results in 


dx _ m gz(t) 
(5.3) dt sim—1)+1P. P 


which yields upon integration 


m ‘ : 
(5.4) x= [s@m—1) +1, P J avoae 


Assuming the over-all injection velocity as constant, the equation for the 
characteristic lines can be expressed as 


mt 
(5.5) x= is@n—1) +172 X const. 


[s(m 
where the ‘‘const.”’ depends only on the porosity and the injection velocity 
(and is thus independent of the viscosity ratio). The characteristics, thus, are 
straight lines; they are shown in Fig. 6. 

Assuming that the displacement in the porous medium starts with a sharp 
front, one can now calculate its spread due to fingering. Since the characteristics 
are straight lines in the x—¢ plane, the speed c(s) with which each saturation 
value proceeds, is given by the slope of that line: 


(5.6) c(s) = const. lena 


where the ‘‘const.’’ is the same as that in eq. (5.5). The speed for each satura- 
tion value is independent of time; thus the distance x(s) to which it travels 


is given by c(s)t: 
m 
(5.7) x(s) = const. is(m—1)+1P° t. 


This is the connection between x and s for the time ¢ which is sought after. 
A sketch of it is given in Fig. 7. 

The manner in which time is involved in eq. (5.7) shows that no stabiliza- 
tion is possible in the present model of fingering. The basic form of the satura- 
tion curve (Fig. 7) simply grows ad infinitum by similarity stretching. 
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It may be of interest to rewrite the above formulas in terms of pore volumes 
(é) and pore lengths (@). Thus, considering a system of length L, one may 
introduce the following dimensionless variables: 


(5.8) ff aeat /(PL) 


(5.9) x/L. 
Then one has from (5.4) 


m .. 
(5.10) = [s@m—1)+1F é. 


The break-through value of ¢ is obtained if s is set equal to zero, and Z equal 
to 1; hence 


(5. 11) Leena _ 1/m. 


Similarly, the value of f for complete recovery is obtained if s is set equal to 
1 and & equal to 1; hence 


(5.12) tcomplete recovery = ™. 


It is obvious from eqs. (5.11) and (5.12) that the above development can 
only have meaning when m 2 1. 


NOMENCLATURE 


(s) = velocity of saturation movement (cm sec“). 
permeability (cm?). 

relative fovrnend ee dimensionless). 
length a a system (cm) 

viscosity ratio (dimensionless). 

pressure (g cm™ sec™?). 

porosity (fractional; dimensionless). 

filtration velocity (cm sec“). 

fractional flow (dimensionless). 

fictitious saturation (fractional; dimensionless). 
= time (sec). 

= velocity (cm sec™). 

= linear co-ordinate (cm). 

= viscosity (g cm™ sec™). 


Ss 


| 
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EQUIVALENCE BETWEEN CONTINUOUS AND 
DISCRETE RADIATING ARRAYS! 


A. KSIENSKI 


ABSTRACT 


The radiation patterns produced by continuous excitation distributions and 
discrete arrays are compared and the conditions are derived under which one 
type of source may be substituted for the other with negligible errors. It is shown 
that the aperture lengths in both cases should be the same but the element 
spacing is dependent on the type of pattern desired. Examples are computed to 
demonstrate these relations for both directive patterns and shaped beams. 


INTRODUCTION 


The pattern synthesis problem received considerable attention in the past, 
and a large amount of information is presently available. Many of the results 
are given for continuous distributions; in practice, however, it is hard to 
control amplitude and phase independently across a continuous aperture. 
Usually, then, the continuous excitation distribution is used as the envelope 
for the excitations of a discrete array. It is the purpose of the present paper 
to examine the relationship between continuous and discrete arrays and 
estimate the errors involved in substituting one for the other.* It is shown 
that the current practice of using the continuous distribution as an envelope 
is justified in most cases. [t is also shown that the parameter to be adjusted 
for different classes of patterns is the element spacing and that the aperture 
width should remain the same if the error tolerance is to be preserved. Examples 
are computed to show these effects for patterns of greatest interest, namely 
directive beams and cosecant squared patterns. 


THE EQUIVALENCE THEOREM 


Let F(u) represent the pattern produced by a linear aperture with a con- 
tinous distribution (u = mw sin 6, where @ is the observation angle). Let P(u) 
represent the pattern produced by a linear array of point sources. It will be 
shown that under certain conditions P(u) = F(u). If the array represents a 
Fourier series it can also be shown that the feeding coefficients A, are pro- 
portional to the continuous distribution at corresponding points. 

Let F(u) be a function prescribed for —x < u < mw under the following 
conditions: 

No. 1. F(u) is of bounded variation (which is satisfied by almost all prac- 

tical patterns) ; 

No. 2. F(u) = 0 for |u| > mr. 


‘Manuscript received September 6, 1960. 

Contribution from the Aerospace Engineering Division, Hughes Aircraft Company, Culver 
City, California. 

*Several of the mathematical results derived in this paper may be obtained from Shannon's 
sampling theorem. It was felt, however, preferable to derive all the results in a consistent form 
in antenna terminology rather than to transcribe formulas derived in another field and attempt 
to justify their applicability, particularly in view of the simplicity of the mathematics involved. 
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The second condition is closely approached by a large number of directive 
patterns. Let the continuous excitation distribution which produces F(u) be 
given by g(p) where ? is the spatial, or position variable, measured in half 
wavelengths (p = x/(A/2)). Then (Silver 1949) 


(1) F(u) = J elode™p, 


(2) (6) = sc J Flue ™du. 


However, if F(u) = 0 for u > a, then 


(3) eo) = 5 - Fe Mau. 


Now consider an array of discrete sources which are represented by an ortho- 
gonal set of functions, whose orthogonality is defined over the range from 
—7x to m. For example, a set of 2NV+1 point sources, where consecutive sources 
are displaced by a distance \/2 from each other, the resultant pattern of 
such an array will be given by 


(4) P(u) = . Af” 


where u = m sin 0, and where A, are the amplitude and phase of the excitation 
coefficients. In order to approximate F(u) by means of P(u) we have to 
determine appropriately the coefficients A, of P(u). The coefficients are given 
by 

1 . —jnu 
(5) A = xf F(uje "du. 


2a 


The right-hand side of equation (5) is identical with the one of equation (3) 
except for m replacing p. Hence, 


(6) A, = g(n). 


Thus it is seen that the feeding coefficients A, for the discrete array are 
obtained from the continuous distribution g(p) by inserting the proper integer. 
Hence for any pattern F(u) obeying condition No. 2, its Fourier transform 
forms the envelope of its Fourier coefficients. 

No mention has been made yet of the size of the aperture, or the total 
number of sources. Assume the continuous aperture producing F(u) is M 
wavelengths wide, namely 


g(p) = 0 for |p| > M. 


Actually g(p) can be only approximately zero for p> M, but M being 
arbitrary, it may be so chosen as to make g(p) as close to zero as desired 
for p > M. (Otherwise the integral in equation (1) would not converge.) We 
will show now that the necessary discrete array is also M wavelengths wide. 
Using equation (6) 
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A, = g(n) 
and from the above condition on g(p) it follows the 
A, = 0 for |n| > M 


or all feeding coefficients beyond the Mth will be zero. This requires that the 
array have 2M@-+1 sources, or be M wavelengths long, i.e. the same as the 
continuous aperture. 

We shall presently prove that the pattern produced by the discrete array 
would be exactly F(u) if the condition that g(p) = 0 for |p| > M were satis- 
fied. Given a complete orthonormal set of functions {f,(u)} defined over the 
interval —x < u < m, and given F(u) on the same interval and subject to 
condition No. 1 the summation >> ynf,(u) converges to F(u) where y, are given 
by (Courant and Hilbert 1953) 


(7) n= f F(u)fn(u)du. 
The set {e+*“/./2r} where n— © is a complete orthonormal set, hence 
(8) F(u) = 3 Ane™ 
where 
(9) An = Yn/V 2. 


However, all A, are zero for |n| > M, hence 


(10) F(u) = Y Age™= Plu), 


which proves that the pattern has been exactly reproduced by the discrete 
array. Note that F(u) was specified only over the interval |u| < 7, and the 
discrete array of 2M-+1 elements will reproduce F(u) exactly over that 
interval. The equality between F(u) and P(u) breaks down, of course, out- 
side the visible range that is, for |u| > a, due to the essential difference between 
the Fourier series and the Fourier transform. It might be well to point out 
that one cannot have both F(u) and g(p) vanish identically outside a finite 
range of the variables, hence the assumption that g(p) = 0 for |p| > M 
implies that F(u) does not identically vanish for |u| > x although its magni- 
tude may be as small as desired. This means that if one expects to approximate 
very closely a desired pattern by means of a finite aperture one has to refrain 
from specifying F(u) to vanish outside the visible range. In this respect the 
discrete array is somewhat less restricted since in the latter case one has to 
request only that g(p) vanish at integral values of p. 


Error Evaluation 

If condition No. 2 is satisfied only approximately an error will result in 
the excitation distribution that will be proportional to that part of the integral 
which was neglected in the Fourier expansion. More often the practical 
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approximation problem will be of the following nature. The continuous 
aperture which is finite in extent will usually not be sufficient to contain the 
distribution g(p) required for the prescribed pattern. In other words, the actual 
span of the aperture is only a portion of that required by an exact synthesis. 
For example, the required length is 1 wavelengths where M is very large and 
the available aperture is only N wavelengths wide, where N < M. This results 
in an error in the desired pattern F(u). Let us denote this resulting pattern 
by F.(u). The discrete array of the same span will likewise produce an 
approximate pattern P,(u). However, the error produced by the discrete 
array will not necessarily be larger than that of the continuous aperture. If 
we denote the mean square error caused by the continuous aperture as EF, 
and that of the discrete array as Eq we shall obtain for the discrete array, 


7 N 


Pe 2) ‘ a oa 2 ee ‘2. ox jnu ; 
(11) E, = i f |F(u) —P,.(u)|"du = a: F(u) do Ane du 


1 7 N 


an J. | F(u)|"du 2. A, 


then 
—N M 


(12) E, = a lanl + 2 |A,|’. 


For the continuous aperture 


aT 7 N 2 
(13) E, = < a |F(u) — F,(u)|"du = - J | F(u)— J. g(pye™dp| du 


27 


7 ° N 
as if i 3 |F(u)|*du+ f : \F(u)l'duh —2 7 lg(p)|*dp. 


Using Parseval’s relations in equation (13) we obtain 


(14) F=f leorltae— flee) lap 


where we have assumed that F,(u) is negligibly small for |u| > z, i.e., that 
its contribution to the total error is small. This assumption is justified in the 
majority of patterns, particularly the directive ones. From equation (14) we 
can see that the error in the case of the continuous distribution is proportional 
to the part of the area of |g(p)|? which was left out due to the finite extent 
of the aperture. In the discrete distribution the error is proportional to the 
sum of the |A,|? which were left out due to the finite extent of the discrete 
array. Clearly the two errors are very close to each other. In the case that 
F.(u) is not negligible in the range |u| > a additional terms have to be 
included in the error expression 


(ia) Fe= J le@)iap— J le lap—se J Feu) fd 


| 

—M N 
+ fru) lau 
ro c e\(uU ° 
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The last two terms reduce the error of the continuous aperture in comparison 
to the discrete one. This reduction is actually a correction term since the first 
two terms in equation (14a) represent the total r.m.s. deviation of F,(u) from 
F(u) for all values of u over the infinite range, while E, is the r.m.s. error over 
the finite limits +2. Thus it is clear that in the case where the continuous 
aperture produces a pattern which significantly extends into the invisible range 
the discrete array of half wavelength spaced elements will produce a poorer 
approximation. This case is discussed in more detail in the following sections 
where the appropriate inter-element spacing is derived in order to approximate 
more closely the continuous aperture. 

We now consider an example of the case where condition No. 2 is only 
approximately satisfied (Fig. 1). Then we compare in Figs. 2 and 3 the patterns 


CONTINUOUS DISTRIBUTION ——— 
DISCRETE DISTRIBUTION —— — 
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Fic. 1. Patterns produced by a Taylor excitation distribution. Comparison between a 
continuous distribution 5\ long and an array of 11 elements spaced half wavelength apart. 


of uniformly illuminated continuous apertures and corresponding discrete 
arrays which illustrate the case when both only approximate a desired pattern. 


Examples 
Taylor Pattern 
A Taylor pattern (Taylor 1953) 





CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


CONTINUOUS DISTRIBUTION —— 
DISCRETE DISTRIBUTION — — —~ 
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Fic. 2, Patterns produced by a uniform illumination. Comparison between a continuous 
distribution 5d long and an array of 11 elements spaced half wavelength apart. 


sin r/(u"—B*) es 

RRS -a6r. — =— . = — 
rV/(u'-—B°) r - 2 

satisfies condition No. 2 approximately, and the patterns produced by the 

continuous distribution and the discrete array should be approximately the 

same. The continuous distribution is given by 


(16) g(b) = Jol jBV (xe —p’)]; p= ax/a 


where B = 0.7386 and the aperture width = 5A. Now compare this with a 
discrete array which consists of 11 source points, at \/2 separation, so the 
over-all distance is 5A. The excitation coefficients are 


mend Lole-))) 


As is seen from Fig. 1 the patterns are quite similar to each other. The 
main beam of the discrete array is somewhat narrower and the sidelobes are 
shifted correspondingly. The physical reason for this behavior is that the 
“space factor’ effect of a point source is equivalent to a half wavelength 
aperture centered at the point source. Hence an array consisting of NV point 


(15) F(u) = 
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CONTINUOUS DISTRIBUTION 
DISCRETE CISTRIBUTION ——— — 


RELATIVE ELECTRIC FIELD 
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Fic. 3. Patterns produced by a uniform illumination. Comparison between a continuous 
distribution 6A long and an array of 12 elements spaced half wavelength apart, i.e., 54A long. 


sources spaced \/2 apart may be considered to be equivalent to an array 
N/2 wavelengths long rather than (N—1)/2 which is the physical length of 
the array. The main difference lies in the element factor of a half wavelength 
continuous aperture which multiplies the pattern of the continuous aperture 
and thus progressively reduces the sidelobes. The mathematical explanation 
can be obtained directly from the error formulas, equations (12) and (14). 
Thus imagine an aperture distribution which vanishes at a distance of (N+3) 
half wavelengths from the origin. From equation (12) it is clear that an 
array of 2N+1 sources would produce the exact pattern since all terms after 
the Nth would vanish, thus the array would be N wavelengths long. On the 
other hand the continuous distribution g(p) would have to extend to (N+) 
half wavelengths on each side of the origin to produce the exact pattern. This 
would result in an aperture whose length would be (N+3)A, or one half 
wavelength longer than the discrete array. From these considerations it is 
clear that this effect may result in a difference of up to one wavelength. (This 
limit depends on the orthogonality interval for which the Fourier expansion 
is defined.) 





342 CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


Uniform Illumination 

Another rather classical case is the uniform illumination. Let us examine 
and compare the pattern produced by the continuous distribution versus the 
corresponding discrete array. A uniformly illuminated aperture of over-all 
dimension a produces the field pattern, 


= pele ne -l<u<l 


(18) F.(u) 
while for a discrete array of N half wavelength spaced elements the pattern is 


sin[(NV/2) ru] 
sin[(4/2)u] ° 


The reason that the subscript (e) was attached to the above pattern expres- 
sions is due to the fact that the uniform illumination is the result of attempting 
to produce a delta function whose Fourier transform is a constant. However, 
due to the finiteness of the aperture the excitation is cut off at a finite distance 
from the origin. Thus both the continuous and the discrete distribution are 
finite approximations to the required distribution. Consider a continuous 
aperture 5A long and uniformly illuminated and an array of the same length 
with its element spaced \/2 apart thus consisting of 11 elements. From Fig. 2 
one can see that the resulting patterns are quite similar except for the effects 
that occurred in the previous example, that is, the main beam is narrower in 
the discrete case and the sidelobes shifted accordingly. As an additional 
example consider a 6A continuous aperture and compare it to 53A discrete 
array, both uniformly illuminated. In the present case (see Fig. 3) the two 
patterns are almost identical and the only difference is in the sidelobe ratio. 
The difference between the patterns shown in the above example is a good 
indication of what one might expect in substituting a discrete array for a 
continuous aperture. It might be worth noting that the good agreement 
between the patterns is due to the fact that all of the above patterns approxi- 
mate quite closely condition No. 2. 


(19) P.(u) =B -l<u<l 


Energy Consideration 

Observe that condition No. 2 refers to the invisible range (|u| > x) of 
u = msin 6. The behavior of F(u) in that range is related to the reactive 
energy stored in the vicinity of the antenna (Woodward and Lawson 1948). 
Thus the physical significance of condition No. 2 is the requirement that the 
antenna be predominantly resistive, namely, that it should be an efficient 
radiator. It may be interesting to note that whenever Fourier coefficients 
are computed, the integration is carried out only over the visible range 
(|sin 6| < 1). Thus if the pattern of the continuous aperture has a large value 
of F(u) in the region |u/z| = |sin | > 1 it would not carry over into the 
discrete array. This implies that an array representing a Fourier set of functions 
is of low reactive nature, and, for example, free from supergain effects.* 

*The reactive energy would, of course, be infinite if one were to consider infinitesimal 


elements, since the pattern would repeat indefinitely in the invisible range. The finite size, 
however, of physical elements insures that the pattern vanish for u>1. 
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THE EQUIVALENCE FOR PATTERNS NOT SUBJECT TO CONDITION NO. 2 


For most cases the assumption of F(u) = 0 for u > = is quite realistic, 
that is, F(u) is either zero or close to zero. This condition implies directive 
beams, for which the pattern is desired to vanish outside a small interval of 
u covering the main beam. We shall presently extend the treatment to non- 
directive patterns which require a non-vanishing signal for the end of the 
visible range, namely a certain amount of “end fire’’ energy. For such patterns 
F(u) usually vanishes outside an interval +-ow where o > 1. We will show 
that if the function F(u) is non-zero for intervals larger than z, the spacing 
of the elements in the array necessary for good pattern agreement decreases 
to values smaller than \/2. This indicates mathematically that the array is 
susceptible to increased amounts of reactive energy. The same is clear from 
physical considerations due to the increased coupling between elements caused 
by their proximity to each other. As the interval for which F(u) ¥ 0 tends 
to infinity the spacing tends to zero, producing a continuous distribution. The 
converse is also true, namely, if F(u) = 0 outside an interval of u smaller 
than 7 the spacing between the elements may be larger. The faster F(u) 
vanishes the larger can the spacing between the elements become. This would 
permit the use of fewer elements in the antenna without increasing the error 
of approximation. This last effect, however, is not very useful, since for 
spacing much larger than /2 the pattern will repeat in the visible range, 
producing a main beam or a portion of it in the sidelobe region. 

Consider F(u) = 0 for u > om where o > 1. Using the same notation as 


before 


s 3 ®oo , sami 7 1 f 4 ies 
(20) ge) =p J Rwemau =f Fwemdu 


In order to obtain relations similar to the previous ones between the discrete 
and continuous excitations, we consider a set of functions w, defined to be 
orthogonal over the interval form —omr to +o7. Then 



























f w, . w,du = 0 for k #1 


—oTr 






and 





oT 


(21) A aw, .w,du = | fork =1 


—OoTr 












where A is a constant independent of & and /. 
Let us consider a linear array of radiators and determine how the above 
requirement can be met. The pattern of such an array is given by 







(22) P(u) = > A,e™ 





where u = msin 6, and a is the distance in half wavelengths between con- 
secutive radiators. To approximate a prescribed pattern F(u), let 
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F(u) =~ >> A,e™ 


a ee F(uje7 "du 
2a —7 /a : 


Note that if we let a = 1/o we obtain 


(24) A, = J. F(uje "du, 


Qna 


which is the desired result, and w, = e”“/”, The above relation means that 
the distance between consecutive sources has to decrease from \/2 by the 
same factor as the range of F(u) * 0 extends beyond +z to permit the 
synthesis of the pattern in the invisible range as well as in the visible 
portion. As F(u) extends more into the invisible range, the spacing of the 
sources decreases producing more coupling, and acquiring more of the dis- 
advantages of continuous distributions and making the physical construction 
of the array almost as hard, if not harder, than the continuous distribution. 


Example 

To give an example of a pattern that benefits by spacing the array elements 
at distances somewhat smaller than \/2 we now consider a pattern which is 
required to have a non-vanishing signal at one end of the visible range (u = 7) 
and zero signal at the other (« = —7). The desired pattern is specified as 


Fu) == =2 ae ee 


(25) u 


F(u) =0 —@r<qgu < 2/13 x. 


Considering an array of elements spaced 4/2 apart we obtain for the array 
expression 


N 
P(u) = >> A,e™, 
—N 


which is periodic for intervals of 27. This implies that 
Pix) = P(-s) 


and thus the conditions given by equation (25) cannot be simultaneously 
satisfied at u = +. The resulting pattern, using a Woodward synthesis, would 
compromise the two conditions and the signal would fall somewhere between 
unity and zero. This is shown in Fig. 4 by the broken line. The present 
difficulty can be easily removed by extending the synthesis interval into the 
invisible range. Thus if 15 elements are placed within the same space of 6A 
as was formerly occupied by the 13 elements, the resulting inter-element 
spacing will be [(13/15)\/2] and the synthesis range or the interval of periodi- 
city will extend to « = +15/13 7. The result is shown by the solid line of 
Fig. 4. Pi(u) is given by 
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PATTERNS 
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Fic. 4. A cosecant squared pattern obtained by a Woodward synthesis for a discrete 
array 6\ long. Comparison between the pattern resulting from 4/2 spacing and that resulting 


from an inter-element spacing .434). 
; sin| 13 (. ae) 
7 a “Aa 
115 2 13 


(26) P,(u) = dX a 
“(2 3(4-2)] 
716° 9 13 


and P,(r) = 0.7961; Pi(—2) = 0.4221. The value of P:(u) for u = x is thus 
within approximately the same error as for the rest of the range, while the 
sidelobes in the zero signal range are below the same peak value. The modified 
pattern has somewhat smaller ripples near u = 7; the reason being that the 
first zero signal point has been removed further from the non-zero portion of 
the pattern. The invisible range can thus serve as a region into which some 
of the difficulties of the visible range synthesis may be carried and settled. 
In the present example the change in reactive energy is negligible, however, 
a useful objective has been achieved by using a set of elements spaced some- 
what less than A/2 apart. 

In the above example the continuous aperture excitation required to produce 
the pattern described by the solid line would have a non-vanishing pattern 
outside the visible range. The failure to notice it would produce a distorted 
pattern in the visible range if one would use the continuous excitation distri- 
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bution as the envelope for the discrete array. To emphasize the importance 
of checking the total pattern produced by a continuous distribution consider 


the following pattern, 


(27) F(u) = 0 —-r<qu<d0 
= 1, O<ucr 


which may be synthesized equally well by either a continuous aperture or a 
discrete array of elements spaced \/2 apart. Computing the continuous 
excitations we obtain 

pirP/2) sin (px/2) 


es, 
(28) (6) = 5c [oman = a 
The discrete array results in 


(29) A, = elf" sin (mm /2) 
ne 
This approximation may be quite satisfactory except at u = 0 and u= 7 
due to Gibb’s phenomenon. 
The vicinity of u = 7 may be easily adjusted in the continuous aperture 
by adjusting the excitation to extend the pattern into the invisible range simi- 
larly to the previous example. This would involve a change of the excitation 


which would be 


1 ar : Pe 6 ‘ 
(30) g(p) = +f e Ip des om e' p/2) a / ) 


where 1 < a < 2. If this excitation were used as an envelope for the discrete 
array we would obtain 


= phar) Sin awn 


(31 ) An rn 


This excitation would result in the following pattern for a 4/2 spaced array 


sin aan Ta 
(32) Pe): = ae nl a 


which would extend the main beam into the range of angles between —7 
and zero: 


P(u) = —t Qu < (a—2)r 
(a-2)r cu <0 
Ocucr 


In other words a beam of width am is produced whose center is scanned to 
u = anr/2. Thus a completely unacceptable pattern is produced unless one 
uses element spacing of \/2a. This obviously does not mean that an acceptable 
pattern cannot be obtained by means of an array of elements spaced \/2 
apart. It only means that in this case the continuous distribution cannot be 
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used as the envelope for the discrete array. The pattern has to be directly 
synthesized for the given array since the inter-element spacing was pre- 
specified. 


EQUIVALENCE THEOREM FOR TWO-DIMENSIONAL DISTRIBUTIONS 


The one-dimensional derivation can be generalized to include two-dimen- 
sional planar apertures. We can represent the far field of a two-dimensional 
aperture by means of a transform (Silver 1949) 


(33) g(kz,Ry) es . £ F(é,n)e2 ="! dtdn 


and 


1 2. 3 
(34) Fen) = gia Jf elbebye ated, 
where we shall measure the variable £ and 7 in half wavelengths as before, 


Me 


and 
k, = rsin@cos @¢ 


(35) — 
ky = wsin @sin ¢. 


Here g(k:,k,) is the desired pattern usually specified for the visible range 

alone, namely for —2 < k,,k, < x. The conditions imposed on g(k;,k,y) are: 
condition No. 1, g(kz,k,) is of bounded variation as a function of k, and k,; 
condition No. 2, g(kz,ky) = 0 for |k,,ky| > 7. 


The procedure for proving the equivalence between continuous distributions 
and discrete arrays is the same as in the one-dimensional case,* and will, 


therefore, not be repeated. 
The error formula is likewise similar and is given by 


; 1 © © . 
(36) Ey = zh ff le(bcky)|"dbcdky— I [Ara 
where 


1 . 7 —jmkr.—jnky 
Ame ™ G53 if i g(ksky eo "dhe edhey. 


CONCLUSIONS 
We have shown that a discrete array can always be substituted for a 


continuous distribution in a plane, without changing the pattern significantly. 
This holds for both directive and shaped beams. The over-all dimensions of 


*For details of obtaining a two-dimensional Fourier series representing a planar array, 
which is required for the derivation, see (Ksienski 1960). 
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the array are the same as that of the continuous distribution, and the excitation 
coefficients of the array are proportional to the continuous excitations. The 
spacing of the sources in the array are determined from the behavior of the 
pattern outside the visible range. Increasing the density of sources per unit 
length beyond the one indicated by the above considerations does not improve 
the resulting pattern. Certain minor differences in the pattern may occur, 
depending on the orientation of the elements in the array. Thus, for a given 
antenna dimension, a set of transverse slots would produce the effect of a 
somewhat longer and less directive array than a line source. This is evidenced 
by a somewhat narrower beam and higher sidelobes produced by the discrete 
array. For a set of longitudinal slots, however, these differences would be 
virtually eliminated. 
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APPENDIX 
ERRORS IN APPROXIMATING A CONTINUOUS FINITE APERTURE 
BY A SMALLER APERTURE, DISCRETE OR CONTINUOUS* 


We are interested in evaluating the mean square error that results when the 
pattern, F(u), produced by a continuous linear aperture of finite length 2M 
is approximated by that produced by a continuous linear aperture of length 
2N < 2M. This is in contrast to the assumption made in the main text that 
the aperture size for obtaining F() exactly is infinite. The present assumption 
results in F(w) extending over an infinite range, the error being evaluated only 
over the visible range. The aperture distribution of the smaller aperture is 
identical with that of the larger but is merely chopped off at +N. Denote the 
aperture distribution as g(p) where p is a normalized aperture variable and 
call the pattern from the shortened aperture F,(u). Then the mean square 
pattern error in the visible range is given by 


(1a) By = 5 J (Fw) — Feu) IP (u) — Fe(u)] du. 


*Derived by Dr. A. T. Villeneuve, Hughes Aircraft Company, Culver City, California. 
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If the patterns are expressed in terms of their aperture distributions the 
error may be put into the following form 


- - «/~\ sin (p—f) er 
Qa) Ee= ff etoyer ee) SP=O* apar 


~ «/.\ sin (p—f) 
+ ff sore) SPEDE spar 


eT * * sin (p—{) 
— fi J, eve) +e"oreey] LRD* par. 


This may be written in a form as follows 


- «/.\ sin(p—f) xr 
(3a) E.= Jf J eee PSO aa 


where A is given by four squares of side |MM— | in the (f,p) plane located 
at the four corners and inside the large square whose corner points are given 
by (M, M; —M, M; —M, —M; M, —M). 

If we attempt to approximate F(u), the pattern from a continuous aperture 
of normalized length 2M, with a discrete array of elements separated by 
ad/2 the resulting pattern will be periodic with period 27/a and will coincide 
with F(u) in the interval —2/a < u < w/a. Call the periodic pattern F;(x). 
Then 


(4a) Fi(u) = D> ane™™ 


The source distribution which gives rise to this may be written formally as 
1 a —jpu 1 7 j(neu—pu) 

(5a) gi(o)=5-} Filuje ™ du = 5 > ane du 
= —o —ca 


=} a,5(p—an). 


This represents a discrete aperture distribution of point sources where in view 
of equation (4a) and the conditions preceding it, 


m/f 


(6a) Qn = -" (uje "™ du = =f F(uje”™™ du. 


2a —T/a 
Due to the finite available aperture the number of discrete sources are limited 
to 2N + 1 where N is the largest integer in M/a. Therefore there will be an 
error in the resultant pattern. If the mean square error between F(z) and the 
resultant pattern is computed one gets 


am M M 2 sin(p—¢) 
(7a) Ee= ff etover) SPO apar 


*sina(m—n)r 


N 
+d Tate a(m—n)r 
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DIPOLE-FIELD TYPE MAGNETIC DISTURBANCES 
AND AURORAL ACTIVITIES! 


B. K. BHATTACHARYYA? 


ABSTRACT 


The characteristics of the magnetic field components at Agincourt have been 
calculated for a current system produced by an electric dipole located in the region 
of auroral activity near Ottawa. It is noted that, irrespective of the orientation of 
the dipole, the horizontal magnetic field component rotates in the clockwise and 
anticlockwise senses for motion of the dipole towards the east and the west 
respectively, when the dipole is situated in the north half of the sky as seen from 
the observing station. 

Next, the magnetograms obtained at Agincourt have been studied at those 
times of the night when auroral activity was recorded in the all-sky camera 
photographs at Springhill near Ottawa. It is noted that the horizontal magnetic 
field describes a loop during a particular phase of auroral activity because of its 
gradual growth and decay. The distributions of clockwise and anticlockwise 
rotations with respect to local time are found to be very similar in many respects 
to those of auroral motions to the east and west respectively. The sense of rotation 
of the loop is predominantly anticlockwise in the early part of the night and clock- 
wise in the late hours of the night. 

It is found that eastward and westward orientations of the dipole are the most 
probable ones. The direction of movement and the initial location of the pre- 
dominant auroral form in the sky are found to tally well with those of the dipole 
deduced from a study of the magnetograms. 

It seems that there is a time sequence relationship between successive phases of 
auroral activity and changes of characteristics of the loops described by the hori- 
zontal magnetic field vector. The area of a loop and the maximum magnitude of 
the field vector in the loop appear to be related to the brightness and horizontal 
extent of the auroral forms. 


INTRODUCTION 


From the early days of research on auroral and geomagnetic phenomena, it 
has been known that the onset of auroral activity is very closely associated in 
time with the magnetic perturbations recorded on the surface of the earth. 
This has led many to investigate whether the current system responsible for 
magnetic perturbations is spatially coincident with the visible auroral forms 
in the sky or not. ‘Polar elementary storms’ which are practically always 
accompanied by aurora, were shown by Birkeland (Harang 1951) to be attribu- 
table to the field of an infinite linear current or to that of a fairly narrow sheet 
of current flowing in the upper atmosphere. Since then attempts have been 
made to correlate the location of the linear current and that of the auroral 
arc observable in the sky. 

The lines of the field of force of an infinitely long current conductor are 
represented by a system of circles around the conductor and so the magnetic 
disturbance vector will be tangent to one of these circles passing through the 
place of observation. This fact can be utilized to determine the angle of eleva- 
tion of the infinite linear current from the horizon by knowing the direction of 
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the disturbance vector. Stagg and Paton (1939) have calculated for several 
cases the elevation y of the auroral arc and the elevation 6 of the line current. 
They found that the frequency of y and 6 lying on the same side of the magnetic 
vertical plane was five times that on the opposite side. The difference between 
6 and y was less than 10 degrees in 35 per cent of the cases and less than 20 
degrees in 60 per cent. A similar conclusion was reached by Heppner (1954) 
from his analysis of magnetic and auroral records at College, Alaska. 

A study on the same lines was undertaken at Ottawa with the help of auroral 
all-sky camera photographs taken at Springhill (geographic 45.2° N., 75.5° W.; 
geomagnetic 56.5° N., 6.9° W.) and the magnetograms taken at Agincourt 
(geographic 43.8° N., 79.3° W.; geomagnetic 55° N., 11.6° W.). Since Spring- 
hill is about 230 miles northeast from Agincourt and the latitude separation 
is only 1.4°, it can reasonably be assumed that the effect of current associated 
with auroral forms observed at Springhill will be felt at Agincourt. 

First the all-sky camera photographs are studied for the three months of 
March, April, and June, 1958, which cover the Spring equinoxial period when 
auroral activity was at its peak in that year. Some nights on which the sky 
condition was perfectly clear and the elevations of auroral arcs could be 
reasonably accurately measured, are chosen. Then the magnetic records on 
those nights are scaled at intervals of 15 minutes to obtain the values of the 
horizontal component H (positive northward), the declination D (positive 
eastward), and the vertical component Z (positive vertically downward) of the 
earth’s magnetic field. With the help of the median values Hy, Dg, and Z, of 
the 5-quiet-day values in the respective month, the deviations AH, AD, and 
AZ are calculated. AH and AD are then combined together trigonometrically 
to determine the disturbance horizontal magnetic field vector, Hy. Hg and AZ 
are used to find the disturbance field vector. 

Though the preliminary calculations indicate the validity of the conclusion 
reached by Stagg and Paton on many occasions, there are, nevertheless, 
periods ranging from half an hour to several hours on every night studied 
when the perpendicular to the disturbance field vector points consistently deep 
into the south whereas auroral activity is confined within 45° of the northern 
horizon. A typical example is shown in Fig. 1 for a 5-hour period on the night 
of April 16, 1958. During this period arcs and rayed arcs, confined within 30° 
of the northern horizon, constituted the auroral activity. But from the Hg 
and AZ plots in Fig. 1, it is evident that the infinite ‘line’ current, if responsible 
for the perturbed magnetic field, was always in the southern part of the sky. 
If we take the effect of earth current into account and assume that twice the 
vertical component AZ and two-thirds of the horizontal component AH are of 
external origin, the position of the ‘line’ current is shifted still further to the 
south. The line current model of the magnetic disturbances does not therefore 
associate itself spatially with the visual auroral arc during a considerable 
length of time on every night at Springhill near Ottawa. Even when the 
locations of the line current and auroral arc are in the same part of the sky, 
they differ, in most of the cases, by quite an appreciable angle in elevation as 
was noted by Stagg and Paton. 
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Fic. 1. Magnetic field disturbance vectors for a period of 5 hours on April 16, 1958. 


In addition to what has been mentioned, there is another factor to be noted. 
Bands, draperies, glow, and such other forms which are limited in horizontal 
extent from the observational point of view, are sometimes the main con- 
stituents of auroral activity. During these periods, it does not seem reasonable 
to assume an infinite linear current or many such currents as representing the 
perturbed magnetic field. 

Furthermore, if a series of vectors of the horizontal magnetic field com- 
ponent are plotted on a polar diagram, the position of the end of the vector 
rotates clockwise or anticlockwise and forms a loop during any active phase 
of auroral activity. The description of a loop by the horizontal component 
cannot be due only to the growth and decay of the intensity of the current 
system but must also result from some kind of movement or rotation of the 
current system as a whole. An infinite line current obviously cannot explain 
all the features of this characteristic of the magnetic field disturbance. 

In recent years Nagata and Fukushima (1952) and Fukushima (1953) have 
shown that polar magnetic storms are composed of a number of elementary 
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disturbances which take place intermittently or successively with duration 
from a few tens of minutes to about 2 hours. It has been suggested that each 
elementary disturbance is caused by the production of an electromotive force 
in a rather small area of high ionization produced in the upper atmosphere 
over the auroral zone. The equivalent current system corresponding to each 
elementary disturbance is generally much simpler, being approximated by the 
current system produced by an electric dipole at the center of the highly 
ionized area. 

The present study is undertaken to determine as far as practicable the 
suitability of representing the cause of the magnetic field disturbances by a 
current system due to an electric dipole spatially located in the region of 
brightest auroral activity. An explanation of the loop characteristics of the 
horizontal magnetic field disturbances will be attempted by assuming the 
movement of the electric dipole itself, i.e., of the small region of high con- 
ductivity. Since in the case of severe magnetic disturbances many elementary 
disturbances play their part successively one after another and sometimes 
simultaneously, it is very difficult to separate their individual effects and to 
analyze them one by one. In this paper, therefore, only those disturbances, 
which are not very complicated and are classified under small or medium range, 
have been studied. We have considered only the case when a single electric 
dipole is present in the sky. 


MAGNETIC FIELD COMPONENTS OF AN OVERHEAD CURRENT SYSTEM DUE 
TO AN ELECTRIC DIPOLE 
When an electric dipole directed to A = 0 is situated at the pole of a spherical 
co-ordinate system (0, A) of a conducting spherical surface, the current 
function J of the resulting electric current is given by (Fukushima 1953) 


Toe sin O sin A 


(1) J(0,4) = 2a ° 1—cos 9 





where a denotes the radius of the spherical surface of the ionosphere, Jo the 
moment of the electric dipole, and go the electric conductivity of the spherical 
surface outside the dipole. This expression of the current function is valid 
everywhere except in the neighborhood of the region of high conductivity 
containing the dipole. Expressing (1) in terms of spherical harmonics, we have 


P,,1(cos @)sin A 





Iooo — 2n+1 

a een n(n+1) 
where P,; is the associated Legendre function of the mth degree and of the 
first order. When the current function is expressed in terms of spherical 
functions, e.g., J = > J,(0,A), the magnetic potential W at the point 
(r,8,A) is given by (Chapman and Bartles 1940) 


(2) J(9,A) = 


n+1 . 
(3) W = ir ee ys, a(Z ) 
(4) = 1204, > at: \’p Py,1(cos 8)sin A 
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where 7 is inside the sphere of radius a and the current and current function 
are expressed in electromagnetic units. 
The generating function of the Legendre polynomial is given by 


1 — i 
where P,,(x) is the Legendre polynomial of the mth degree. 
With the help of the relation 


(6) Pale) = (Inat)? Px) 


we easily obtain 
n—1 ek. 1/2 1 
(7) > h*"P,.i(x) = (1—x’) (Dh FR)? * 


n=l 


Integrating both sides of (7) with respect to 4 from 0 to h, we have 


(h—x) 
(8) . E Pale) = Ga aan. eal 


With the help of (4) and (8) we can express the magnetic potential W in the 
following form: 


a Loo 1 f (h—x) | 
(9) = ‘9. Al yt G thx th)! Zs 


where x = cos 9 and h = r/a. 

We now introduce a new spherical co-ordinate system in which the axis 
passes through the geomagnetic poles. Let (0,A) be the co-latitude and longi- 
tude of the point of observation in the new system, and let (00,A0) be corre- 
sponding co-ordinates for the position of the dipole. If the dipole be oriented 
geomagnetically eastward, the mutual relations between the two co- ordinate 
systems are 

sin 8 cos A = sin 6 sin (A — Ao), 
(10) sin 8 sin A = cos @ sin 0) — sin 6 cos @) cos (A — Xo), 
cos 8 = cos 8 cos 0) + sin 6 sin 8) cos (A — Apo). 
If the dipole be oriented along a geomagnetic meridian towards the south, 
we have 
sin 8 cos A = sin 6 cos @) cos (A — Ao) — cos 8 sin Bo, 
(11) sin 6 sin A = sin 6 sin (A — Apo), 
cos 8 = cos 6 cos 6) + sin 8 sin 0) cos (A — Ao). 


With the aid of (10) and (11) we can express (9) in terms of geomagnetic 
co-ordinates and then determine the northward, eastward, and vertically 
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downward components of the magnetic field in the usual way. The expressions 
of these components for a dipole oriented in the geomagnetically east direction 
are: 

Northward component: 

ow = 4a Ioao = A 


(12) X= oe aes — AJ — sin 0 cos 6-008 @ 


sin 0) cos(A— nt . C+{sin @ sin 60-+cos 6 cos 8 cos(A—Ao)} . | 
Eastward component: 


sin 0) sin(A—Ao)C 


a ase _ OW 4 Jed _ sin A 


rsin@d\ 7° 2a sin @ 
+008 6 sin(x—)B | 


Vertically downward component: 


_ 2. a9 fate __hAsin Asin 9@_ 
(14) a- Da -2. 2 


Similarly, for the dipole oriented to the south, we have 
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r ° 2a sin 0 


\ —sin 6 cos 09+ cos @ sin 09 cos(A— } 


. C—cos 6 sin(A— Xo) . B| 


«a Jas sin A 
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(16) 


- 4m fe _hsin Asin @ 
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where 
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In equations (12)-(17), the trigonometrical functions of A and © are used 
for the sake of writing simplicity. The first two equations in both cases of 
dipole orientation are for geomagnetic north and east components respectively. 
They can be easily manipulated to calculate the expressions for the geographic 
north and east components. Then these latter components are combined to 
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yield the horizontal field component H and the declination D. In Fig. 2 are 
plotted the expressions for (ar/2xIoa0)H (solid lines) and (ar/24Iya9)Z (dotted 
lines) as a function of various locations and two different orientations (geo- 
magnetic east and south respectively) of the electric dipole taking the 


N 


GEOMAGNETIC LONGITUDE ——> 


Fic. 2. A plot of the theoretical magnetic field components at Agincourt as a function of 
various locations of the electric dipole and for two different dipole orientations marked E 
(geomagnetic east) and S (geomagnetic south) respectively. At the top of the figure N, E, S, 
and W stand for the four cardinal directions, the geographic north, east, south, and west 
respectively. The solid lines indicate (ar/2Ioa0)H and the dotted lines (ar/2xIoa0)Z. 


observing station at Agincourt (@ = 35° and A = 348.4° E.). The dipole is 
assumed to be at a height of 100 km from the surface of the earth. 

From the expressions of the magnetic field components for two perpendicular 
orientations of the dipole, it is quite easy to obtain their values for any other 
dipole orientation, whatsoever. For example, if the dipole be oriented west- 
ward, the directions of H and Z will be opposite to those in the case of east 
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orientation. For a dipole pointing at an angle with the geomagnetic east, the 
values in cases of south and east orientations will have to be combined 
vectorially to obtain the field components. 

A study of Fig. 2 reveals some very interesting facts. The horizontal com- 
ponent is found to rotate either clockwise or anticlockwise as the dipole moves. 
The direction of the rotation depends upon the location of the dipole and the 
direction to which it is moving. In the following table are summarized the 
various possible cases for specific sense of rotation of H. 


Direction of 
Sense of rotation movement of ‘ c Orientation of 
of H the dipole Location of the dipole the dipole 


(i) Clockwise East North of the parallel of latitude All possible 
through Agincourt 
West South of the parallel of latitude All possible 
through Agincourt 
North West of the geomagnetic All possible 
meridian through Agincourt 
South East of the geomagnetic All possible 
meridian through Agincourt 


(ii) Anticlockwise West North of the parallel of latitude All possible 
through Agincourt 
East South of the parallel of latitude All possible 
through Agincourt 
South West of the geomagnetic All possible 
meridian through Agincourt 
North East of the geomagnetic All possible 
meridian through Agincourt 


Aurora is usually observed in the northern sky at Springhill except for very 
disturbed days. Since Agincourt is only 1.5° south of Springhill, these aurorae 
will also be on the same part of the sky at Agincourt. Auroral motion is pre- 
dominantly in the geomagnetic east-west direction, though a small north- 
south component of motion is always present (Kim and Currie 1958; Bhatta- 
charyya 1960). There are some cases when the motion is mainly in the north— 
south direction. If the region of high conductivity in the upper atmosphere 
which contains the electric dipole be spatially coincident with the brightest 
and most active auroral form seen in the sky, the dipole will also move in most 
cases along the geomagnetic east or west direction. Consequently the horizontal 
field component due to current system of the dipole will rotate either in 
clockwise or anticlockwise sense and will eventually form a loop because of the 
gradual growth and decay of the intensity of the current system. 

It is, indeed, difficult and susceptible to error to predict the motion of the 
dipole from the sense of rotation of the horizontal field component. Neverthe- 
less, by noting the sense of rotation and the mean direction of H and the 
variation of the vertical field component it is often possible to deduce the 
direction of motion, orientation and location of the dipole and to associate it 
with the important auroral form in the sky as recorded in the all-sky camera 
photograph at that time. 
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CHARACTERISTICS OF THE LOOP FORMED BY THE HORIZONTAL FIELD 
VECTOR DIAGRAM 


In the previous section we have already discussed the reason why the vector 
diagram of the horizontal magnetic field component due to an electric dipole 
whose moment increases up to a peak value and then decreases in a gradual 
fashion with time, will form a loop. The sense of rotation of the horizontal 
disturbance vector will depend upon the location and, most important of all, 
the direction of motion of the dipole. 

In this section we shall study the characteristics of these loops formed on 
those nights when aurora was recorded on the all-sky camera at Springhill 
during the months of March, April, and June, 1958. The magnetograms used 
for this purpose were obtained from Agincourt and only small and moderate 
intensity disturbances were analyzed. 

Generally, the magnetic activity follows the auroral activity. The latter 
undergoes a cycle of changes in its form, illumination, spread, and motion. 
Any particular period during which one specific form builds up to its full 
grandeur and then decays gradually due to diffusion, de-excitation, and a 
host of other factors, may be called a phase of auroral activity. It should be 
noted that a phase does not in general consist of one specific form but several 
different forms and that it may not be isolated in time, i.e., before one phase 
has died down, one or more other phases may start. Due to this overlapping 
of phases, the loop described by the horizontal disturbing force will, in general, 
not be complete. So it is normally expected that before a particular vector 
diagram has reached even less than the half-way mark in its decay, a new 
phase will start. This tallies with what is actually observed. 

Since the duration of a phase, and hence that of a loop, ranges from a few 
tens of minutes to 2 hours, we have scaled the magnetograms normally at 
every 5- or 10-minute interval. The vector diagrams for the horizontal com- 
ponent of the magnetic field are plotted in Figs. 3 and 4 for selected periods of 
high auroral activity. It is evident from the plots that the individual loops can 
be identified, though they are not necessarily isolated from other loops. The 
period of a particular loop can also be determined from the plots without 
difficulty. 

Next, the time at which the horizontal field vector in a particular loop 
reaches its maximum, is noted and thus the distributions of clockwise and 
anticlockwise rotations of loops with respect to time are obtained. The dis- 
tributions are shown in Fig. 5. The number of anticlockwise rotations is quite 
high in the early evening and assumes peak values during the intervals 0100- 
0200 and 0400-0500 U.T. After local midnight, the number begins to decrease 
gradually. In the case of clockwise rotations, the number begins to increase 
from about 0100 U.T., reaches its peak value in the interval 0700-0800 U.T., 
and then begins to decrease. After 0700 U.T., the number of clockwise rotations 
far exceeds that of anticlockwise rotations. 

It is clear from Fig. 5 that the sense of rotation of the loop is predominantly 
anticlockwise in the early part of the night and clockwise in the late hours of 
the night. The reversal of the sense of rotation from anticlockwise to clockwise 
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Fic. 4. Vector diagrams of the horizontal magnetic field component during specific 
intervals on April 18 and June 15, 1958. The field strength in gammas has been plotted against 
azimuth measured eastward from geographic north. The time for each vector is given in U.T. 
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Fic. 5. Distribution of clockwise and anticlockwise rotations of the horizontal magnetic 
field component with respect to the time of the night. 

Fic. 6. Distribution of the mean angular velocity of the horizontal magnetic field com- 
ponent in degrees per 10 minutes with respect to the time of the night. 
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seems to be a systematic process, occurring at a time very close to local mid- 
night. The distributions of clockwise and anticlockwise rotations with respect 
to local time follow quite closely in characteristics the distributions of auroral 
motions to the east and the west respectively (Kim and Currie 1958; Bhatta- 
charyya 1960), though the times of the peaks of the distributions do not 
coincide exactly. This is to be expected because, as mentioned earlier, the 
east-west component of auroral motion is always accompanied by a normally 
small north-south component. There are a few cases when the latter is quite 
predominant. We have not made any attempt to separate the east-west 
component of motion from the north-south one. This is the reason why the 
effect of the north-south component of auroral motion will be somewhat 
reflected in Fig. 5. It seems that most of the cases of clockwise rotations before 
0400 U.T. are due to southward progress of auroral forms located to the east 
of the geomagnetic meridian passing through Agincourt. 

From the disturbing force vectors at times just before and after the time of 
occurrence of the maximum horizontal field vector in a particular loop, the 
mean angular velocity of the rotation of the vector for all loops in a particular 
hour is calculated. The mean angular velocity is plotted in Fig. 6 as a function 
of universal time. It is found that the angular velocity is maximum within 
approximately the interval 0300-0400 U.T. for anticlockwise rotation and 
within 0400-0500 U.T. for clockwise rotation. We note (Bhattacharyya 1960) 
that the approximate time interval is 0200-0300 U.T. for the occurrence of 
maximum speeds in the case of auroral motion to the west and 0300-0400 U.T. 
for that to the east. It seems that the angular velocity of the field vector de- 
scribing the loop lags somewhat behind the corresponding auroral motion either 
in the eastward or westward direction. After midnight, the angular velocities 
in both the curves of Fig. 6 are gradually decreasing. 


AURORAL MOTIONS AND MOTIONS OF THE CURRENT SYSTEM 


Let us now attempt to deduce the direction of motion, orientation and loca- 
tion of the dipole by comparing the sense of rotation of the horizontal field 
vector and the nature of variation of the corresponding vertical field component 
observed at Agincourt with the relevant theoretical curves in Fig. 2. 

The distribution of orientation thus obtained is as follows: 50 cases of east- 
ward orientation, 32 westward, 8 northward, and 18 southward. It seems that 
eastward and westward orientations are the most probable ones. Strangely, 
no relationship between the orientation and direction of motion of the dipole 
has been noticed. As an example, an eastward-oriented dipole has been found 
to move in practically all directions. 

It is encouraging to note that at most of the times when auroral motion is 
clearly detectable, the direction of movement and initial location of the pre- 
dominant auroral form in the sky tally well with those of the dipole. Three 
examples are described in the following table. In Figs. 3, 4, and 7 are shown 
the horizontal field vector diagrams and the all-sky camera photographs for 
three nights in April, 1958. 








PLATE [| 


Fic. 7. All-sky camera photographs taken at Springhill on the nights of April 17, 18, and 
30, 1958. The serial day number appears at the upper left of each photograph and the U.T., in 
hours and minutes, at upper right. 
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Characteristics Possible dipole Characteristics 





; Sense of 
Date and time rotation of in the variation orientation and of auroral 
interval the loop of Z location motion 
(1) April 17, 1958 Anticlockwise; Z always East-oriented An auroral band 


0450-0520 U.T. H always 
negative and 
minimum at 
0510 U.T. 


(2) April 18, 1958 Clockwise; H 
0410-0450 U.T. _ negative and 
minimum at 


0420 U.T. 


(3) April 30, 1958 Clockwise; H 
0500-0645 U.T. negative and 
minimum at 


0610 U.T. 


positive; maxi- dipole; moving 
mum at westward; 
0510 U.T. initial location 
N.E. of 
Agincourt 
Z always South-oriented 


dipole; moving 


positive and 
south; location 


maximum at 


0430 U.T. east of 
Agincourt 
Z always West-oriented 


dipole; moving 
east; location 
N.W. of 
Agincourt 


negative; 
minimum at 
0610 U.T. 


moving west- 
ward from N E, 
horizon to 
western 

horizon 


Auroral curtain 


mainly located 
in the east; 
moving 
towards south 


Auroral patches 


mainly located 
in the N.W.; 
moving 
towards east 


A study of cases like those given in the above table suggests that there is a 
good possibility of the highly conducting region being coincident with the 
moving auroral forms. 


AURORAL ACTIVITY AND CHARACTERISTICS OF THE HORIZONTAL FIELD 
VECTOR DIAGRAM 


Previously we have discussed to some extent the various phases of auroral 
activity which range in period from a few minutes to half an hour or more. The 
durations of successive loops formed by the horizontal field vector also vary 
within the same range. Some sort of relation may then be reasonably expected 
between the different phases of auroral activity and the corresponding loops 
of the horizontal magnetic field. With this in mind, auroral activity is studied 
on different nights from all-sky camera records and compared with the diagrams 
of the horizontal field vector. It is, indeed, noticed that (a) most of the loops 
are associated with particular phases of auroral activity; (0) the width and 
extent of a loop increase with the intensity of illumination, horizontal extent, 
and movement of the particular forms taking part in the corresponding phase 
of auroral activity; (c) the overlapping loops are observed when the decay of 
one phase and the growth of the succeeding phase are not separated by an 
interval of time. 

Different phases of auroral activity on two nights are described below and 
the corresponding horizontal field vector diagrams are shown in Fig. 4. 

The following examples definitely indicate a time sequence relationship 
between auroral activity and characteristics of the loop described by the 
horizontal field vector. The area and maximum magnitude of the loop seem 
to have a connection with the brightness and horizontal extent of the auroral 


forms. 
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Characteristics Possible dipole Characteristics 





Sense of 
Date and time rotation of in the variation orientation and of auroral 
interval the loop of Z location motion 
(1) April 17, 1958 Anticlockwise; 2Z always East-oriented An auroral band 
0450-0520 U.T. 4H always positive; maxi- dipole; moving moving west- 
negative and mum at westward; ward from N E. 
minimum at 0510 U.T. initial location horizon to 
0510 U.T. N.E. of western 
Agincourt horizon 
(2) April 18, 1958 Clockwise; H Z always South-oriented Auroral curtain 


0410-0450 U.T. negative and 
minimum at 


0420 U.T. 


(3) April 30, 1958 Clockwise; H 
0500-0645 U.T. negative and 
minimum at 
0610 U.T. 


positive and 
maximum at 
0430 U.T. 


Z always 
negative; 
minimum at 
0610 U.T. 


dipole; moving 
south; location 
east of 
Agincourt 


West-oriented 
dipole; moving 
east; location 


N.W. of 
Agincourt 


mainly located 
in the east; 
moving 
towards south 


Auroral patches 
mainly located 
in the N.W.; 
moving 
towards east 


A study of cases like those given in the above table suggests that there is a 
good possibility of the highly conducting region being coincident with the 
moving auroral forms. 


AURORAL ACTIVITY AND CHARACTERISTICS OF THE HORIZONTAL FIELD 
VECTOR DIAGRAM 

Previously we have discussed to some extent the various phases of auroral 
activity which range in period from a few minutes to half an hour or more. The 
durations of successive loops formed by the horizontal field vector also vary 
within the same range. Some sort of relation may then be reasonably expected 
between the different phases of auroral activity and the corresponding loops 
of the horizontal magnetic field. With this in mind, auroral activity is studied 
on different nights from all-sky camera records and compared with the diagrams 
of the horizontal field vector. It is, indeed, noticed that (@) most of the loops 
are associated with particular phases of auroral activity; (b) the width and 
extent of a loop increase with the intensity of illumination, horizontal extent, 
and movement of the particular forms taking part in the corresponding phase 
of auroral activity; (c) the overlapping loops are observed when the decay of 
one phase and the growth of the succeeding phase are not separated by an 
interval of time. 

Different phases of auroral activity on two nights are described below and 
the corresponding horizontal field vector diagrams are shown in Fig. 4. 

The following examples definitely indicate a time sequence relationship 
between auroral activity and characteristics of the loop described by the 
horizontal field vector. The area and maximum magnitude of the loop seem 
to have a connection with the brightness and horizontal extent of the auroral 


forms. 





CANADIAN JOURNAL OF PHYSICS. VOL. 39, 1961 


Auroral activity 





(a) April 18, 1958 


From 0107 rays and arcs are visible in the northern sky. 


At 0112 rays become very bright. At 0121 they have 
been transformed into glow and arc which begin to de- 
crease in intensity rapidly after 0126. They regain their 
intensity of illumination after 0136 and break up into 
separate glows and bands at 0140. At 0145 homogeneous 
glow covers the whole northern sky right up to zenith, 
spreads southward, and thus becomes diffused. 


At 0153 a band appears in the N.E. part of the sky, 


growing in intensity and moving westward. At 0156 
several bands and glows cover the northern sky up to 
zenith. At 0212 glows, bands, rays, and arcs fill up the 
sky up to 35° S. of zenith. This activity begins to recede 
towards the north at 0216 and diminishes in intensity. 


At 0221 new bands and arcs appear and extend up to 


15° S. of zenith. The spread is maximum —45° south of 
zenith at 0237. 


At 0330 several auroral forms are visible up to 32° south 


of zenith, become very bright, and stay up to 0348. 
From 0349 separate bright bands and glows appear in 
the north and are later transformed into glow and 
homogeneous arcs confined to the northern sky practi- 
cally up to zenith. 


At 0421 a very intense band is formed in the northeast, 


moving towards southwest, becoming brightest at 0423, 
and breaking up at 0427. This activity subsides down 
to a great extent at 0439. 


At 0446 two quite bright bands are growing in intensity, 


one in the eastern horizon and the other in the north- 
west sector of the sky. The former is moving towards the 
west and spreading in the north-south direction at 0450 
whereas the latter is moving towards the south and also 
spreading in the north-south direction. At 0454 the 
former is very bright and extends from the northern 
horizon to 60° south of zenith through the zenith. It 
gains in intensity and spreads with time. At 0504 it is 
decreasing in brightness and has broken up into two 
parts. At 0514 quite bright auroral forms are again 
observed in the N.W., N.E., and S.W. sectors of the sky, 
moving north at 0516. At 0530 there is no trace of this 
auroral activity. 


Cloudy intervals. 


At 0552 auroral activity is again building up in the whole 


northern sector of the sky up to zenith. Auroral forms 
are found to move to the north and are reduced to a 
narrow arc in the northern horizon at 0600. At 0608 a 
few bright patches and bands appear in the north. They 
are transformed later into a bright glow which stays up 
to 0626. 


(b) June 15, 1958 


Up to 0530 patches of cloud prevented observation. 
At 0540 auroral forms such as arcs, glows, and bands 


are observable in the sky up to 55°S. of zenith and 
quite active up to 0550. 


Characteristics of loops 


Anticlockwise rotation of a loop 
from 0050 to 0140. The hori- 
zontal field vector is maxi- 
mum at 0130. 


Anticlockwise rotation of a loop 
from 0150 to 0230. The hori- 
zontal field vector is maxi- 
mum at 0210. 


Clockwise rotation of a loop 
from 0330 to 0410. The hori- 
zontal field vector is maxi- 
mum at 0340. 


Clockwise rotation of a loop 
from 0410 to 0440. The hori- 
zontal field vector is maxi- 
mum at 0420. Note the high 
magnitude of the field vector 
at maximum. 


Anticlockwise rotation of a loop 
from 0440 to 0530. The hori- 
zontal field vector is maxi- 
mum at 0510. 


Anticlockwise rotation of a loop 
from 0540 to 0620. The hori- 
zontal field vector is maxi- 
mum at 0610. 


Clockwise rotation of a loop 
from 0535 to 0550. The hori- 
zontal field vector is maxi- 
mum at 0545. 
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Auroral activity Characteristics of loops 


At 0554 auroral activity is concentrated mainly near the Clockwise rotation of a loop 
zenith and consists of a very bright band and a faint from 0555 to 0615. The hori- 
glow to the north. At 0556 the band has moved slightly zontal field vector is maxi- 
to the north and at 0558 has broken up into two seg- mum at 0610. 
ments, one near eastern horizon and the other near 
western horizon of the sky. At 0600 the segment in S.E. 
sector of the sky becomes prominent, continues to move 
toward N.W. through zenith, and begins to decay in 
intensity. At 0608 another band appears near zenith in 
N.W. sector, moving to S.E., expanding, and getting 
diffused. This exists up to 0616. 


At 0618 arcs and glows are observed in the north. At 0620 Anticlockwise rotation of a loop 
a part of an arc is transformed into a bright band ex- from 0615 to 0630. The hori- 
tending from about 60° above the northern horizon to zontal field vector is maxi- 
zenith. It is getting diffused from 0622 and vanishes at mum at 0625. 

0630 leaving only diffuse glow and arc in the north. 


From 0632 another band is growing in intensity near the Anticlockwise rotation of a loop 
zenith. It proceeds gradually towards the south and from 0635 to 0650. The hori- 
extends up to 38°S. of zenith at 0642. After 0642 it zontal field vector is maxi- 
begins to break up and becomes diffused. At 0650 there mum at 0645. 
is no trace of the band. 

Cloudy intervals. 





SUMMARY AND DISCUSSION 


In this paper are presented the results of an attempt to represent the 
magnetic field perturbations associated with aurorae by a dipole field. It has 
been shown that the clockwise and anticlockwise senses of rotation of the 
horizontal magnetic field vector can be explained by assuming the region of 
high conductivity containing the dipole to be coincident in space and moving 
along with the predominant auroral form in the sky. 

The electric dipoles are found to be oriented in all directions, though they 
seem to prefer the eastward or westward directions to the others. These orienta- 
tions do not appear to be related to the directions of velocities and so it does 
not seem likely that they are produced by a dynamo action of the highly 
ionized area, as suggested by Nagata and Fukushima (1952) and Fukushima 
(1953). Furthermore, for the dynamo action to be valid we have to invoke a 
close connection between the wind motion of the upper atmosphere and auroral 
motion. But normally auroral forms move with a very high velocity with 
respect to the ionospheric wind motion and this does not tend to justify any 
general relation between them. 

The assumption of a single dipole may not hold good on many occasions 
when there are several auroral forms simultaneously active in the sky. This 
is the reason why severe magnetic disturbances have not been analyzed in our 
study. In the case of disturbances of small or moderate intensity, the presence 
of a single dipole seems to be able to explain most of our observations. If there 
are several dipoles simultaneously active but moving in the same direction, the 
horizontal magnetic field will rotate in the same sense as in the case of a single 
dipole. Hence the distributions of clockwise and anticlockwise rotations with 
universal time (Fig. 5) need not be modified. However, in cases when these 
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conditions are not satisfied, it is not at all simple to predict the locations and 
directions of movement of a number of dipoles by comparing observational 
curves with theoretical ones. 

It has been noticed in this study that magnetic perturbations are normally 
associated with a display of rays, draperies, and other active forms. They are 
very rarely detected when there is only a quiet homogeneous arc in the vicinity 
of the northern horizon. These observations may be explained by the discharge 
theory of aurora (Reid 1958) which states that at the base of the E-layer of the 
ionosphere, there will be very little extra ionization associated with homo- 
geneous arcs but large amounts of ionization with the onset of an active display 
of rays, draperies, pulsating and flaming aurorae. 

Recently, Weaver and Skinner (1960) have derived the current systems 
induced by dynamo action in a sheet of partially ionized gas surrounding an 
elliptical region of greater ionization density and have applied the results to 
conditions existing in the ionosphere in the presence of a long, homogeneous 
auroral arc. It has been shown that in the case of a long arc, a line current whose 
magnitude is linearly dependent on the height-integrated Cowling conductivity 
within the arc, can account for the observed magnitude of the disturbed mag- 
netic field. A line current model has also been used by Bless et al. (1959) to 
explain the magnitude of the bay-type disturbances associated with aurora. 
Though a line current may, in some cases, be quite a reasonable model just 
for estimating the magnitude of the magnetic field disturbances, it cannot 
explain, as noted in this study, the clockwise or anticlockwise rotation of the 
magnetic field vector. Also in cases of auroral forms, limited in horizontal 
extent, the assumption of a line current model seems hardly to be justified. 
It appears, therefore, that a dipole model is, in general, preferable to a line 
current model for a correlation study of the different characteristics of the 
magnetic field disturbances and of the successive phases of accompanying 


aurora. 
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TYPICAL FAILURES IN PULSED FIELD MAGNET COILS* 
RICHARD STEVENSON 


A pulsed field magnet capable of producing fields in the megagauss range 
has been used here for some time. It appears worth while to describe, in the 
light of our experience, some of the factors which affect its successful operation. 

The magnet itself is patterned after the type built at the M.I.T. Lincoln 
Laboratory (Foner and Kolm 1956). Energy is supplied to the magnet from 
a 2000-uf bank of 3000-v capacitors. The capacitors are charged by a separate 
power supply, and discharged into the magnet coil by means of an electrical 
spark. The maximum energy available from the capacitor bank is 10,000 
joules; a field of 500,000 gauss can be produced with about 2600 joules in a 
solenoid (described below) of a 3/16 in. inner diameter. 

The basic element of the magnet is a flat helical solenoid machined from a 
rod of beryllium copper with a parting tool, and heat-treated. The standard 
designs used here are described in Table I. In the assembly of the coil, adjacent 

















TABLE I 
Inner diameter Outer diameter Thickness of Approximate useful 
(in.) (in.) No. of turns turn maximum field 
3/8 1} 10 .060 300,000 
1/4 1 10 .060 550,000 
3/16 1 10 .060 750,000 





turns are separated from each other by washers of shock resistant impregnated 
teflon or of Fiberglas impregnated with teflon. Our experience is that the 
reinforced teflon is superior to mica insulation for this purpose, and far more 
convenient. 

The magnetic forces which act on the coils are of two types. There is a 
radial force tending to increase the inner diameter, which is a maximum at the 
center. Also there is an axial force which tends to warp inwards the innermost 
portion of each turn (Furth and Waniek 1956). These forces are resisted by 
the mechanical strength of the conductor. Photographs of a coil which has 
failed under these circumstances are given in Figs. 1 and 2. This coil was 
subjected to approximately a megagauss; at the center, the inner diameter 
was doubled, and the turns were fractured at each end. The failure is, of 
course, explosive, and apparatus placed in the solenoid is destroyed or severely 
damaged. 


*Supported by the Office of Naval Research and in part by the Defence Research Board 
(Canada). 
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Fic. 1. Notice the bulging of the sides which is characteristic of a failure at very high fields. 

Fic. 2. This is the same coil as in Fig. 1. The inner diameter has doubled, and the originally 
flat turns have been warped inward. The force of the field has also fractured the coil sym- 
metrically. 


This type of failure always occurs when the magnetic field goes beyond 
700,000 to 800,000 oersteds. Careful heat treating will increase the maximum 
field slightly, but eventual failure is inevitable. 

An entirely different type of failure is characteristic of lower field operation. 
The inner diameter of our magnet coil is covered carefully with a ceramic 
cement, which then is molded with a paper straw. The ceramic-straw assembly 
is allowed to harden for several days. In an experiment the maximum field is 
varied over a wide range. We find that eventually small bulging flaws appear 
on this inner coating. Shortly after the appearance of these flaws the magnet 
will explode, even at rather low fields. 

Figure 3 is a photograph of a coil which has failed in this manner. The metal 
of this coil retained its mechanical strength. (Continued application of the 
high fields appears to anneal the beryllium—copper, even though the heating 
effect is negligible. This in itself appears to set a limit on the lifetime of a coil.) 
The failure seems to be due to sparking, which has cut a slot in the metal of 
about 2 mm deep and $ mm wide running in a line from turn to turn from the 
center to one end. 

It is believed that this type of failure is due to a small flaw in the coil, 
perhaps a tool mark on the inner diameter. When the magnetic field is applied, 
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Fic. 3. This photograph gives a view of a coil which has failed at relatively low fields. The 
failure is probably due to a tool mark on the inner diameter. Careful surface preparation of the 
coil after heat treatment gives an improvement in the lifetime. 


presumably a stress concentration appears at the flaw. This would try to 
relieve itself by propagating a crack radially. This in turn would cause a large 
voltage gradient at the end of the crack, and a spark could occur which would 
damage the insulation and would flaw the adjacent turn. The cracking— 
sparking effect is cumulative and the flaw runs to the end of the coil. 

Careful attention to the machining and surface preparation seems to improve 
the life of the coils. A good coil can withstand several hundred pulses, although 
others will fail almost immediately. The lifetime of a coil is roughly inversely 
proportional to the highest field at which it is used. All coils fail eventually, and 
it is wise to use a fresh coil when starting an experiment. 


Foner, S. and Koi_m, H. H. 1956. Rev. Sci. Instr. 27, 547. 
—— 1956. M.I.T. Lincoln Laboratory Group Report M35, 67. 
Furta, H. P. and Waniek, R. W. 1956. Rev. Sci. Instr. 27, 195. 
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POTENTIAL ENERGY CONSTANTS, ROTATIONAL DISTORTION CONSTANTS, 
AND THERMODYNAMIC PROPERTIES OF STIBINE 


S. SUNDARAM 


INTRODUCTION 


The present investigation on the stibine molecule is a continuation of a 
systematic study by Sundaram, Suszek, and Cleveland (1960), of the hydrides 
of the fifth group elements along with their deuterated and tritiated isotopic 
forms. A preliminary report, on the investigation of the spectrum of SbH;, 
has been made by Smith (1951). It indicated the necessity for accurate high 
dispersion measurements to identify the individual lines and to treat the 
vibration-rotation problem in detail. Haynie and Nielsen (1953) have made 
such dispersion measurements on SbH; and SbD;. They have corrected the 
band centers for anharmonicities so as to estimate with considerable accuracy 
the values of the normal frequencies. In the present investigation, the complete 
quadratic potential energy function for the stibine molecule and its isotopic 
forms has been determined from the vibrational spectral data. 


POTENTIAL ENERGY CONSTANTS 


The details, with regard to evaluating uniquely the six symmetrized force 
constants corresponding to a; and e types of vibrations of pyramidal XY; 
type molecules, have been given by Sundaram e¢ al. (1960). The same pro- 
cedure has been followed here using the harmonic wave numbers of SbH; and 
SbD, given by Haynie and Nielsen (1953). Using the analytical expressions 
relating the symmetry force constants (F) to the constants (f) in the quadratic 
potential energy function, the latter are completely determined. Table I 


TABLE I 
Vibrational spectral data for SbH;, SbD;, and SbT; 








SbH;* SbD;4 SbT;° 
v w v w v w 
ay 1891 1989 1359 1409 1118 1153 
782 796 561 568.5 466 469 
e 1894 1974 1362 1403 1125 1153 
831 845 593 600 486 490 








“Reference Haynie and Nielsen (1953). 
>’Wave numbers obtained in the present investigation using the potential energy con- 
stants given in Table II (see text). 


gives the vibrational wave numbers (w) used in obtaining the potential energy 
constants. v corresponds to the observed wave numbers and w to those cor- 
rected for anharmonicity. The constants are summarized in Table II. The 
vibrational wave numbers given for SbT; in Table I are those calculated in 
the present study with the constants in Table II. The values of v for SbT; 
are obtained using the relation 


wo: = vi ( +x, /0:) 
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TABLE II 
Potential energy constants for the stibine molecule?’ 





Stretching and Bending and Interaction 
interaction between interaction between between bending 
stretchings bendings and stretching 
(md/A) (md A/rad?) (md/rad) 
tn 2.20857 Su 0.73975 fr™ 0.06164 
fr® —0.00398 Snn™ 0.01521 fr’ 0.50895 


“The notation used here is the same as in reference Sundaram et al. (1960). 
‘The number of significant figures is necessary to obtain internal consistency in 
calculations. 


where x; and »; are respectively the anharmonicity factor and the observed 
frequency for the hydride, and »;* and w;* are the quantities for the tritiated 
molecule. It has been shown in the previous paper (Sundaram ef al. 1960) 
that this assumption is quite reasonable. The constants given in Table II 
are transferable to molecules of the form SbXY2 (X and Y = H or D or T) 
belonging to C, symmetry and the results on these will be given elsewhere. 


ROTATIONAL DISTORTION CONSTANTS 


The method of Kivelson and Wilson (1952, 1953) has been applied here 
for the first-order perturbation calculation of the rotational distortion con- 
stants for the molecules SbH2, SbD;, and SbT;. The moments of inertia tensor 
derivatives of the form [Jas‘]o were computed. One needs only [Jzz‘Jo, [Jzz'Jo, 
and [J;2‘]o for such symmetric top molecules. The quantities of the form 


lapys = a [Jas‘lo(F') [Sys ‘To 


are obtained using the symmetrized force constant matrix F. Following Wilson 
(1957), the constants D;, D;x, and Dx are calculated for the three molecules 
and these are given in Table III. There are no experimental values available 
in the literature for a comparison to be made. 


TABLE III 
Rotational distortion constants for SbH3;, SbD;, and SbT3 (Mc/sec) 





Molecule Dy; — Dur Dr 
SbH; 2.4715 4.1586 2.6991 
SbD; 0.6275 1.0572 0.6832 
SbT; 0.2832 0.4759 0.3061 





THERMODYNAMIC PROPERTIES 


From a knowledge of the accurately determined harmonic vibrational wave 
numbers, it has been possible to compute the thermodynamic functions of 
heat content, free energy, entropy, and heat capacity by statistical methods 
for these molecules. The harmonic wave numbers for SbT; are those given 
in Table I and the results for the thermodynamic properties based on these 
should be accurate to within 2%. The calculations of all these quantities have 
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been made for the ideal gaseous state, for 1 atm pressure for a rigid-rotor 
harmonic-oscillator approximation. The values are summarized in Tables IV 
to VI, for various temperatures from 100 to 1000° K, including 273.16° K 
and :298. 
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16° K. 


TABLE IV 
Heat content, free energy, entropy, and heat capacity for SbH;* 





T (°K) (H°—E,°)/T —(F°—Eo°)/T 





100 
200 
273.16 
298.16 
300 


500 
600 
700 
800 
900 
1000 


et ee 


qe 
8. 
8. 
8. 
8. 
400 8.94 
9. 
0. 
Q. 
L. 
Ls 
2 


95 
04 
29 
40 
41 


52 
09 
64 
28 
65 


aL 


38 


51 


53. 
55. 
56. 
57. 
59. 


41 
43. 
46. 
47. 
47. 
49. 
81 


94 
48 
21 
26 
75 


59 
19 
74 
99 


24 


“The values are in cal deg~! mole~! for 1 atm pressure. 


TABLE V 


=. 


46.36 
51.98 
54.77 
55.61 
55.67 
58.69 
61.33 
63.68 
65.83 
68.02 
69.64 
71.35 


C,° 


9. 


11 


15 
16 


95 


7 

8.52 
9. 
9 


46 
83 
85 


20 
12. 
13. 
14. 
15. 


41 
48 
42 
29 
88 


43 


Heat content, free energy, entropy, and heat capacity for SbD;* 

















TCR) (H°—Eo°)/T —(F°—Eo°)/T ss C,° 
100 7.96 40.53 48.49 8.04 
200 8.32 46.12 54.44 9.58 
273.16 8.85 48.79 57.64 10.99 
298.16 9.05 49.57 58.62 11.45 
300 9.06 49.63 58.69 11.47 
400 9.88 52.35 62.23 13.12 
500 10.66 54.64 65.30 14.46 
600 11.39 56.64 68.03 15.52 
700 12.04 58.45 70.49 16.36 
800 12.63 60.10 72.73 17.00 
900 13.13 61.61 74.74 17.48 
13.59 63 .02 76.61 17.87 


1000 


“The values are in cal deg~! mole~! for 1 atm pressure. 


T (°K) (H°=Eo')/T =—(F°—E£0")/T 
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TABLE VI 


Heat content, free energy, entropy, and heat capacity for SbT;* 
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“The values are in cal deg~! mole~! for 1 atm pressure. 
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Similar calculations of the rotational distortion constants and thermo- 
dynamic properties for the partially deuterated and tritiated molecules will 
be reported elsewhere. Determinations of the mean square amplitudes from 
the force fields for all these pyramidal XY; (X = N, P, As, Sb; Y = H, D, T) 


type molecules are in progress. 
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SENSITIVITY OF VG-1A IONIZATION GAUGE 
CALCULATED FROM THE PROBABILITY OF IONIZATION OF GASES 


S. N. GHosH AND B. N. SRIVASTAVA 


INTRODUCTION 

Among the various gauges used for measuring very low pressures, the 
ionization gauge is found to be the most satisfactory in the range of 10% 
to 10-§ mm of mercury. Different types of ionization gauges have been devised 
by different investigators with a wide range of variation in their construction 
and sensitivity. One of these gauges is the type VG-1A which is a modified 
form of the gauge originally devised by Morse and Bowie (1940) and which is 
manufactured by Consolidated Vacuum Corporation, Rochester, U.S.A. 
Young and Dushman (1945) have measured the sensitivity of different 
ionization gauges for inert gases, Hg, He, and Ne and have found that in 
comparison with other ionization gauges the type VG-1 is the most sensitive. 
In this paper, the sensitivity of type VG-1A gauge for a large number of gases 
obtained by Ghosh and Sheridan (1957 and private communication) is com- 
pared with those calculated from the probability of ionization of gases by 
electron bombardments. The gases investigated are inert gases, He, Ne, Cle, 
CO, NO, COs, and hydrocarbons (CH4, C2H4, CoH¢6, and C3Hs). 

SENSITIVITY OF VG-1A IONIZATION GAUGE 
FROM ION GAUGE AND McLEOD GAUGE READINGS 


In column 2 of Table I, the ratios of VG-1A ionization gauge and McLeod 
gauge readings for different gases as observed by Ghosh and Sheridan (1957 
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TABLE I 
Sensitivity of VG-1A ionization gauge from ion gauge and McLeod gauge readings 





Calculated sensitivity 


Observed ratio of (S) of VG-1A gauge for Sensitivity for the gas 
ion gauge to McLeod 1-ma grid current relative to argon 

Gases gauge readings (microampere per micron) (r 

He 6 3.3 0.13 
Ne 3.2 6.3 0.25 
A 0.81-0.77 25.3 1.00 
Kr 0.50 40.0 1.58 
Xe 0.29 68.6 2.72 
Air 1.0 20.0 0.79 
He 1.8 ae 0.44 
Ne 1.0-0.91 21.0 0.83 
Cl, 2.8 7.14 0.28 
CO 1.42 14.1 0.56 
NO Lv 117 0.46 
COz 1.35 14.8 0.58 
CH, 1.76 11.3 0.45 
C.H, 1.20 16.7 0.66 
CoHe 2.3 8.7 0.34 
C;He 2.2 9.1 0.36 





and private communication) are given. (The ionization gauge was calibrated 
against the McLeod gauge for air.) The connection of the ionization gauge is 
of the type (6) given by Dushman (1949), so that the electrons are collected 
by the grid and the ions to the plate. Since a good average value of the sensi- 
tivity for VG-1A ionization gauge below 1 micron is 20 microamperes per 
micron per ma current (data sheet for the ionization gauge tube type VG-1A, 
Consolidated Vacuum Corporation, Rochester, U.S.A.), the sensitivity of 
this gauge for other gases can be obtained by dividing 20 by the values given 
in column 2; the results are shown in column 3. The values of r (sensitivity of 
the gauge for a particular gas relative to argon) are given in column 4. 


SENSITIVITY OF VG-1A IONIZATION GAUGE 
FROM THE PROBABILITY OF IONIZATION OF GASES 

The sensitivity of the VG-1A ionization gauge can also be obtained from the 
probability of ionization of gases by electron bombardments (Reynolds 1931) 
and is given below. 

In Fig. 1 the probability of ionization for different gases as a function of 
electron velocity, obtained by various investigators (Tate and Smith 1932; 
Smith 1930; Compton and Van Voorhis 1925, 1926; Jesse 1925; and Hughes 
and Klein 1924) is given. Tate and Smith (1932) and Smith (1930) have 
obtained data at 0° C, while Compton and Van Voorhis (1925, 1926) derived 
results at 25° C. Other investigators (Jesse 1925 and Hughes and Klein 1924) 
have probably made observations at room temperature. In plotting the above 
probability curves no temperature correction has been made. 

It should be noted that in the experimental determination of the probability 
of ionization of gases per cm path by electron bombardments, the electrons of 
constant velocity traverse the collision chamber only once. On the other hand, 
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Fic. 1. The probability of ionization of gases as a function of electron energy. Full line 
curves are obtained by Tate and Smith (1932) and Smith (1930); dash curves by Compton and 
Van Voorhis (1925, 1926); dot and dash curve by Jesse (1925), and long dash curve by Hughes 
and Klein (1924). 


in an ionization gauge an electron is accelerated in a varying field and makes 
at least one round trip from grid to collector and then back to the grid. 

Assuming the electrodes are concentric cylinders and radial fields linear, the 
field through which the electrons move and the distance traversed by them are 
given in Fig. 2. The dimeasions and normal operating voltages of VG-1A 
ionization gauge as given by the manufacturer are shown below: 
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Fic. 2. Field distribution inside a VG-1A ionization gauge. 


Inner diameter of the collector = 1.375 in 
Diameter of the grid = 9.53 mm 
Filament diameter = 2.1mm 
Filament voltage = 3.0to7.5v 
Grid voltage = +150v 
Grid current = 1 ma 
Collector voltage = —25v 





Most of the ions are produced in the space A’B’ because the electrons after 
passing through the grid acquire sufficient energy to ionize the gas. (The ions 
formed between grid and cathode are collected by the filament and do not 
affect the ionic current.) The voltage at the plane B’ corresponds to the ioniza- 
tion potential of the enclosed gas. As already mentioned, the effective distance 
for ionization of the gas inside the ionization gauge by an electron is at least 
a distance of 2 A’B’ or some multiples of it. 

The amount of ionization of the gas by an electron in the space A’B’ can 
be determined by obtaining the area included by the ionization probability 
curve between ionization potential and maximum accelerating grid voltage 
E, (150 ev), and then by multiplying the result by the distance A’B’. In 
column 2 of Table II, the area of the integrated curve for 1 ma electron (grid) 
current reduced to 1 micron pressure is given. The distance A’B’ for different 
gases is given in column 3. In column 4, the ionic (plate) current 7) for 1-ma 
grid current is obtained by multiplying the values in column 2 and 3. In 
column 5 the average values of 7) obtained by different investigators are 
given. In column 6, the values of 7, obtained by Ghosh and Sheridan (1957 
and private communication) which are the sensitivity (S) of the gauge for the 
different gases are given. The ratio a = 1,/J, which gives the possible number 
of multiples of A’B’ that an electron makes before it is captured by the grid 
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is shown in column 7. Except in the case of CO, NO, and CH, the calculated 
value of a comes out to be about 2, its variation is from 1.8 to 2.7. However, 
since the number of trips which an electron makes between grid and collector 
is an even multiple of A’B’, it follows that in most cases the electron makes 
only one round trip before it is captured by the grid. The calculated values of 
the sensitivity of VG-1A gauge for different gases from the ionization prob- 
ability curve obtained by multiplying the average value of 7) by 2 are given 
in column 8 of Table II. It is to be noted from the last column of Table II that 
Dushman’s (1949, page 350) value for the sensitivity of VG-1 ionization gauge 
is slightly lower. 

The values of r which, as shown before, is the ratio of S for the gas to its value 
for argon were determined previously by Young and Dushman (1945), Lang- 
muir and Jones (1928), and Reynolds (1931). These values and the values 
obtained by Ghosh and Sheridan (1957 and private communication) are given 
in Table III. In columns 6 and 7 of this table, the values of r calculated from 
the ionization probability curve and those given by the manufacturer are 
shown. 


TABLE III 


Comparison of values of r given by different investigators 


r from proba- 





r (R.) bility of ion- 
7({i..J.) r (Y.D.) at 125 ev r (GS.) ization curve 
Gases at 100 ev at 150 ev and250ev_ at 150 ev up to 150 ev r(M.) 
He 0.14 0.13 .10-.13 0.13 .08-.14 0.09 
Ne 0.21 0.21 .12-.24 0.25 .17-.18 — 
A 1.00 1.00 1.00 1.00 1.00 1.00 
Kr — 1.56 —- 1.56 — _— 
Xe — 2.29 — 2.27 —_ — 
He 0.34 0.39 — 0.44 .26-.34 0.28 
Ne 0.90 0.84 .73-.81 0.83 .17-.83 0:7 
CO — — — 0.56 73 0.7 
NO — — 0.56 0.46 .94 — 
CO, — — — 0.58 — 0.74 
CH, “= -- -- 0.45 .68 — 





Note: L.J. Langmuir and Jones (1928). 
Y.D. Young and Dushman (1945). 
R. Reynolds (1931). 
G.S. Ghosh and Sheridan (1957 and private communication). 
M. Manufacturer (Consolidated Vacuum Corporation, Rochester, U.S.A.). 


DISCUSSION 

It is clear from Table III that for He, Ne, A, He, and Nz there isa fair agree- 
ment between the values of 7 calculated from the ionization probability curves 
and those measured by different investigators. The values of 7 obtained by 
Ghosh and Sheridan (1957 and private communication) are roughly the same 
as those obtained by Young and Dushman (1945). 

In the data manual supplied by the manufacturer it is remarked that for He 
and COs, the gauge was influenced by the gettering action (data sheet for the 
ionization gauge tube type VG-1A, Consolidated Vacuum Corporation, 
Rochester, U.S.A.), and hence the sensitivity given by the manufacturer can not 
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be relied on. It is difficult to measure the pressure of these gases with an ioniza- 
tion gauge because these gases react chemically or electrochemically with the 
cathode or other metal elements of the gauge. Also, the McLeod gauge is not 
sensitive enough to give correct readings for compressible gases like CO, NO, 
and CO, at very low pressures. In the case of hydrocarbons, the vapor reacts 
with tungsten filament forming carbides thereby affecting their resistance as 
well as their emissivity. 
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LETTERS TO THE EDITOR 


Under this heading brief reports of important discoveries in physics may be published. These 
reports should not exceed 600 words and, for any issue, should be submitted not later than six weeks 
previous to the first day of the month of issue. No proof will be sent to the authors. 


Frequency Measurement of Standard Frequency Transmissions’: * 


Measurements are made at Ottawa, Canada, using N.R.C. caesium-beam frequency 
resonator as reference standard (with an assumed frequency of 9 192 631 770 c.p.s.). Frequency 
deviations from nominal are quoted in parts per 10". A negative sign indicates that the 
frequency is below nominal. 


GBR, 16 kc/s 
Date, ——————_———- 
November 7-hour 24-hour WWVB, 
1960 average* average 60 ke/s 


N.M. N.M. N.M. 
N.M. N.M. N.M. 
N.M. N.M. N.M. 

N.M. N.M. 
N.M. N.M. N.M. 
N.M. N.M. N.M. 
N N.M. N.M. 
N.M. N.M. N.M. 
N 
N 


SSEE 


2 ie 
a 


v.M. N.M. N.M. 
v.M. N.M. N.M. 
N.M. N.M. N.M. 
— 167 — 166 
— 162 — 164 N.M. 
—165 — 167 
N.M. N.M. N.M. 
—178 N.M. 
—172 —176 
— 166 — 165 
— 163 — 164 
— 161 — 159 
—155 — 159 
— 159 — 160 
— 158 — 157 
—159 — 157 
—153 — 154 
—151 —154 
—154 — 153 
— 148 — 150 
— 148 — 153 
—152 — 150 
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Ke oCoeomOnoar won = 
2S 
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12 
13 
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—_ 
cana 
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Average ¢ — 160 —159 


Midmonthly 
mean of WWV — 149 


Note: N.M. no measurement. 
*Time of observations: 00.00 to 05.30 and 22.30 to 24.00 U.T. 
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